System, method and apparatus for small pulmonary nodule computer aided diagnosis from computed tomography scans

ABSTRACT

The present invention is directed to diagnostic imaging of small pulmonary nodules. There are two main stages in the evaluation of pulmonary nodules from Computed Tomography (CT) scans: detection, in which the locations of possible nodules are identified, and characterization, in which a nodule is represented by measured features that may be used to evaluate the probability that the nodule is cancer. Currently, the most useful prediction feature is growth rate, which requires the comparison of size estimates from two CT scans recorded at different times. The present invention includes methods for detection and feature extraction for size characterization. The invention focuses the analysis of small pulmonary nodules that are less than 1 centimeter in size, but is also suitable for larger nodules as well.

[0001] This application claims the benefit of U.S. ProvisionalApplication No. 60/322,038, filed Sep. 14, 2001, which is incorporatedherein by reference.

BACKGROUND OF THE INVENTION

[0002] The present invention relates to the art of diagnostic imaging ofsmall pulmonary nodules. In particular, the present invention is relatedto analyzing and manipulating computed tomography scans to: segment thelungs, measure lung volume, locate and determine the size of the noduleswithout explicit segmentation, register the nodules using a rigid-bodytransformation, and removing the pleural surface from juxtapleuralnodules in thresholded images.

[0003] Lung cancer is the leading cause of cancer deaths among thepopulation in the United States. Each year there are about 170,000 newlydiagnosed cases of lung cancer and over 150,000 deaths. More people dieof lung cancer than of colon, breast, and prostate cancers combined.Despite the research and improvements in medical treatments related tosurgery, radiation therapy, and chemotherapy, currently the overallsurvival rate of all lung cancer patients is only about 14 percent.Unfortunately the survival rate has remained essentially the same overthe past three decades. The high mortality rate of lung cancer is causedby the fact that more than 80% lung cancer is diagnosed after it hasmetastasized. Patients with early detection of lung cancer followed byproper treatment with surgery and/or combined with radiation andchemotherapy can improve their five-year survival rate from 13 percentto about 41 percent. Given that earlier-stage intervention leads tosubstantially higher rates of survival, it is therefore a major publichealth directive to reduce the mortality of lung cancer throughdetection and intervention of the cancer at earlier and more curablestages.

[0004] The development of the computed tomography (CT) technology andpost-processing algorithms has provided radiologists with a useful toolfor diagnosing lung cancers at early stages. However, current CT systemshave their inherent shortcomings in that the amount of chest CT images(data) that is generated from a single CT examination, which can rangefrom 30 to over 300 slices depending on image resolution along the scanaxial direction, becomes a huge hurdle for the radiologists tointerpret. Accordingly, there is a constant need for the improvement anddevelopment of diagnostic tools for enabling a radiologist to review andinterpret the vast amount of information that is obtained through a CTexamination.

[0005] International Publication No. WO 01/78005 A2 discloses a systemand method for three dimensional image rendering and analysis, and isincorporated herein by reference. The system performs a variety of tasksthat aid a radiologist in interpreting the results of a CT examination.

[0006] One task that radiologists focus on is segmenting the lung regionfrom the image of a single slice obtained from the CT examination. Inthe prior art, some have suggested using a linear discriminant functionand morphological filtering to automatically segment the lungs (S. Hu,E. A. Hoffman, and J. M. Reinhardt, “Automatic Lung Segmentation forAccurate Quantitation of Volumetric X-Ray CT Images,” IEEE Transactionson Medical Imaging, Vol 20, No 6, June 2001, which is incorporatedherein by reference) while others have used mean and median filtering toremove the streaking artifacts due to excessive x-ray quantum noise (J.Hsieh, “Generalized Adaptive Median Filters and their Application inCT,” SPIE, Vol 2299, 1994; J. Hsieh, “Adaptive Trimmed Mean Filter forCT Imaging,” SPIE, Vol 2299, 1994 which is incorporated herein byreference).

[0007] Radiologists also study the location and size of the pulmonarynodules in the CT scan. It is preferred if the radiologist could performthis analysis without the use of explicit segmentation. In some priorwork, the location of a nodule was determined by finding the center ofmass of the nodule through an iterative correlation-based procedure (A.P. Reeves, W. J. Kostis, D. F. Yankelevitz, C. I. Henschke, “Analysis ofSmall Pulmonary Nodules without Explicit Segmentation of CT images,”Radiological Society of North America—2000 Scientific Program, vol. 217,pgs. 243-4, November 2000 which is incorporated herein by reference).The method works for isolated pulmonary nodules, but fails on nodulesattached to the pleural surface.

[0008] Radiologists also estimate a measurement of doubling time of anodule by registering two separate images of the nodule taken at twodifferent times (time-1 and time-2). This analysis requires that thetime-1 and time-2 nodules be registered correctly so that the growth canbe properly measured. Other objects such as vessels and bronchial tubesmust also be registered together. This results in their absence in thedifference image and little effect on the growth measurement.Previously, two nodules were registered by finding the centers of massof the nodules and translating the image accordingly (A. P. Reeves, W.J. Kostis, D. F. Yankelevitz, C. I. Henschke, “Analysis of SmallPulmonary Nodules without Explicit Segmentation of CT images,”Radiological Society of North America—2000 Scientific Program, vol. 217,pgs. 243-4, November 2000 which is incorporated herein by reference).However, this analysis did not guarantee that the two nodules would becorrectly orientated, and that the other objects in the image wouldregistered because these objects might be rotated about the nodule. Somehave registered the nodules by performing a maximization search of themutual information metric over the rigid-body transformation parameters(F. Maes, A. Collignon, D. Vandermeulen, G. Marchal and P. Suetens,“Multimodality image registration by maximization of mutualinformation,” IEEE Transactions on Medical Imaging, vol. 16, no. 2, pgs.187-198, April 1997; Takagi, N.; Kawata, Y.; Nikvi, N.; Morit, K.;Ohmatsu, H.; Kakinuma, R.; Eguchi, K.; Kusumoto, M.; Kaneko, M.;Moriyama, N. “Computerized characterization of contrast enhancementpatterns for classifying pulmonary nodules” Image Processing, 2000.Proceedings. 2000 International Conference on, vol. 1, pgs. 188-191,2000 which are incorporated herein by reference).

[0009] Radiologists also need to remove the pleural surface fromjuxtapleural nodules in CT images. In some prior work, three-dimensionalmorphological filtering and mathematical moments were used to segment ajuxtapleural nodule from pleural surface in a binary image (A. P.Reeves, W. J. Kostis, “Computer-Aided Diagnosis of Small PulmonaryNodules,” Seminars in Ultrasound, CT, and MRI, vol. 21, no. 2, pgs.116-128, April 2000 which is incorporated herein by reference).

SUMMARY OF THE INVENTION

[0010] The present invention is directed to diagnostic imaging of smallpulmonary nodules. There are two main stages in the evaluation ofpulmonary nodules from Computed Tomography (CT) scans: detection, inwhich the locations of possible nodules are identified, andcharacterization, in which a nodule is represented by measured featuresthat may be used to evaluate the probability that the nodule is cancer.Currently, the most useful prediction feature is growth rate, whichrequires the comparison of size estimates from two CT scans recorded atdifferent times. The present invention includes methods for detectionand feature extraction for size characterization. The invention focusesthe analysis of small pulmonary nodules that are less than 1 centimeterin size, but is also suitable for larger nodules as well.

[0011] For the purpose of Computer Aided Diagnosis (CAD), pulmonarynodules are dichotomized into attached nodules and isolated nodulesbased on their location with respect to other solid lung structures.Attached nodules are adjacent to some larger solid structure, such asthe pleural surface. Isolated nodules consist of both well-circumscribednodules and nodules that are larger than all adjacent structures, suchas blood vessels or bronchi. Nodules themselves may be solid, non-solidor part-solid. The analysis of a CT scan for the existence and study ofpulmonary nodules generally entails the following:

[0012] 1. Detection

[0013] (a) Identify the lung regions and main bronchi from thoracic CTimages

[0014] (b) Separate the lungs into two major regions: (1) the lungparenchyma and (2) the lung surface region, including the pleuralsurface and major airways.

[0015] (c) Identify possible locations of isolated nodules in the lungparenchyma region and identify possible locations of attached nodules inthe in the lung surface regions.

[0016] 2. Characterization

[0017] (a) Starting with a single location point within a possiblenodule, identify the nodule region in the CT images. This entailslocating the geometric center of the nodule and approximating its size.

[0018] (b) Given the location and approximate size of a nodule, computecharacteristic features of the nodule, including robust size estimates.

[0019] In connection with the overall methodology of analyzing CT scans,the present invention includes sub-methods for the following:

[0020] 1. The Segmentation of Whole Lung CT Scans into Lung Parenchymaand Lung Surface Regions

[0021] The first preprocessing stage in CT lung image analysis is toobtain the regions of interest from the whole lung scans. The lungregion consists of all tissue found within the pleural surface,including lung parenchyma, vessels, and possibly nodules. Features ofthis approach are the partitioning of the lung into three major regionsand the tailoring of the segmentation algorithm for each region. Inaddition a distinction is made between the central lung parenchyma andthe region of the of the lung parenchyma near to the lung walls. Aproperly segmented lung region greatly reduces the search space of anentire CT scan.

[0022] 2. The Characterization of Lung Air Volume and Inspiration BothFrom the Entire Parenchyma Region and from Single Axial CT Images

[0023] The total volume of the lungs is estimated from the whole lungscans. In addition the change in volume between different scans due tochanges in inspiration is also addressed.

[0024] 3. The Automatic Location and Size Characterization of Nodules

[0025] An algorithm finds the center and approximate size of pulmonarynodules in CT images without the use of explicit segmentation. Thealgorithm has a weighting function that locates the center of mass ofthe nodule. A second weighting function, centered on the location of thefirst weighting function, estimates the size of the nodule. The processis repeated with weighting functions of increasing size until the sizeof the nodule region is reliably determined. The algorithm works forboth isolated nodules and nodules attached to the walls.

[0026] The algorithm is designed for use on images that are resampledfrom high-resolution CT scans (1 mm slice thickness) into an isotropicvoxel space. Such a system could be useful in helping radiologistsanalyze pulmonary nodules. In the CT image browser, the radiologistcould click on a nodule and, using the algorithm, the computer couldautomatically locate the center and size of the nodule. Using thisinformation, the computer then clips out a region of interest for thenodule and performs some nodule characterization analysis and perhapssome 3D visualization.

[0027] 4. The Registration of Nodules from Two Scans

[0028] An algorithm registers two different scans of a nodule region. Athree-dimensional rigid-body transformation of one scan is made tooptimally match the location and orientation (6-degrees of freedom) ofthe second scan. Powell's method is the preferred search strategy.

[0029] The analysis of pulmonary nodules without explicit segmentationestimates a measurement of doubling time by applying a weighted functionon the difference image of the two nodules. The analysis requires thatthe time-1 and time-2 nodules be registered correctly so that the growthcan be properly measured. Other objects such as vessels and bronchialtubes must also be registered together, resulting in their absence inthe difference image and little effect on the growth measurement.

[0030] 5. The Segmentation of Nodules Attached to the Pleural Surface.

[0031] An algorithm segments nodules that are attached to the pleuralsurface from the pleural surface. Starting from a location within anodule, the direction of the pleural surface is determined. A cuttingplane is then iteratively moved towards the surface until the volume ofthe region behind the wall suddenly increases. The increase in volumeindicates that the pleural surface has been reached. The orientation andlocation of the plane is then preferably optimized by a hill climbingprocedure.

[0032] The algorithm removes the pleural surface (and any otherextraneous objects) from a high-resolution image containing ajuxtapleural nodule. The algorithm is not only required to perform agood segmentation, but expected to perform the same segmentation whengiven different scans of the same nodule. This property is necessarybecause the resulting segmented nodule will be used for thecharacterization and analysis of the nodule. If the segmentation is notperformed consistently between images scanned at different times, thenthe doubling time estimation will not be accurate.

[0033] A preferred embodiment of the present invention is a method andapparatus for generating a lung mask for segmenting a lung image fromvoxel data containing the lung image, noise, solid components andsurrounding background information. The apparatus of the invention is amasking unit configured with the teachings of the method of theinvention. The invention can also be practiced on a machine readablemedium. The method includes initially applying a median filter and amean filter to the voxel data to reduce noise. The noise reduced voxeldata is then thresholded to identify solid components and otherstructures. The surrounding background in the noise reduced voxel datais identified and deleted. The connected components of the voxel dataare labeled, and the largest connected components are determined toselect a lung region having a geometric form. The voxel data ismorphological filtered to refine the geometric form of the lung region.

[0034] In a preferred embodiment, the voxel data is associated with aplurality of slices through a patient's lung with the slices beginningat about the patient's shoulders and the application of the medianfilter and the mean filter is performed for the first 25 percent of theplurality of slices only to substantially reduce the computation time.Preferably the median filter has a size of 4×4, and the mean filter hasa size of 1×3. Preferably the voxel data in step is thresholded at agray level of about 500. Preferably the largest connected components areselected to be associated with more than about 1 percent of the voxeldata. Preferably the voxel data is associated with a plurality of slicesfor the morphological filtering with the slices being divided into afirst end region, a middle region, and a second end region. Themorphological filtering is preferably performed with 2D circular filterhaving a first diameter in the first end region and the second endregion while being performed with 2D circular filter having a seconddiameter which is about twice the first diameter in the middle region.

[0035] The present invention also includes a method and apparatus formeasuring lung volume from a segmented lung image. The lung image isobtained from a scan which includes a plurality of slices of voxel datahaving a gray level value and a volume associated therewith. Theapparatus of the invention is a lung volume measuring unit configuredwith the teachings of the method of the invention. The invention canalso be practiced on a machine readable medium. The method includesinitially generating a matrix of entries from the segmented lung image.The matrix includes a plurality of columns and a plurality of rows. Eachof the plurality of columns represent the gray level value and each ofthe plurality of rows represent one of the plurality of slices in thescan with the entries corresponding to the number of times that the graylevel occurs in the corresponding slice. The number of voxels in thesegmented lung image is next determined from the matrix entries. Thenumber of voxels in the segmented lung image is multiplied by the volumeof each voxel.

[0036] The present invention also includes a method and apparatus formeasuring volume of tissue in a segmented lung image. The lung image isobtained from a scan which includes a plurality of slices of voxel datahaving a gray level value and a volume associated therewith. Theapparatus of the invention is a lung tissue measuring unit configuredwith the teachings of the method of the invention. The invention canalso be practiced on a machine readable medium. The method includesinitially generating a matrix of entries from the segmented lung image.The matrix includes a plurality of columns and a plurality of rows. Eachof the plurality of columns represent the gray level value and each ofthe plurality of rows represent one of the plurality of slices in thescan with the entries corresponding to the number of times that the graylevel occurs in the corresponding slice. The sum of tissue of voxels inthe segmented lung image is next determined from the matrix entries. Thesum of tissue of voxels in the segmented lung image is multiplied by thevolume of each voxel. The sum of tissue of voxels is preferablycalculated by summing the product of each matrix entry multiplied by thecorresponding gray level value divided by a gray level value assignedfor tissue.

[0037] The present invention also includes a method and apparatus forfinding the location, P′, and size, r′, of a pulmonary nodule in ahigh-resolution computed tomography (CT) image. The apparatus of theinvention is a nodule finding unit configured with the teachings of themethod of the invention. The invention can also be practiced on amachine readable medium. In the method, a set of initial processingparameters including an initial location, P₁, an initial size, r₁, andtarget value, T, is initially selected. An initial new location P_(i)′of the nodule is computed with a locator template function. Afterincrementally increasing size, r_(i), a new location P_(i)′ of thenodule is computed with the locator template function, and a size metricis computed with a sizing template function. The size, r_(i), isincrementally increased while determining a new location P_(i)′ of thenodule and computing the size metric until the size metric is less thanthe target value. The location, P′, and the size, r′, is then returnedfrom the previous iteration of increasing size, r_(i).

[0038] A preferred method for finding the location, P, and size, r, of apulmonary nodule in a high-resolution computed tomography imageinitially includes windowing the image to ignore bone structures, andselecting a locator template function and a sizing template function. Aset of initial processing parameters including: an initial location, P;size, r; and termination criteria is selected. A search is performed todetermine a maximum response of the locator template function. Aresponse of the sizing template function is determined and compared tothe termination criteria. If the termination criteria has not beensatisfied, the size, r, is incrementally increased, the responsefunctions are determined and compared to the termination criteria. Oncethe termination criteria are satisfied the location, P, and size, r, ofthe nodule are outputted.

[0039] In the preferred method for finding the location, P, and size, r,of a pulmonary nodule in a high-resolution computed tomography image,preferably the image at intensities over about 1000 are clipped towindow the image. Preferably the locator template function is either; aGaussian template function, a Laplacian of the Gaussian templatefunction, or a difference of Gaussians template function. The templatefunction preferably has at least four parameters corresponding to thex-location, y-location, z-location, and radius. The initial location, P,can generally be either calculated from the image or specified by auser. Preferred methods for searching include a hill climbing method andPowell's method.

[0040] The present invention also includes a method and apparatus forregistering 3-d images of a pulmonary nodule from a high-resolutioncomputed tomography (CT) scans. The images include a first image (im₁)obtained at time-1 and a second image (im₂) obtained at time-2, and arein a floating point pixel-format associated with a 6-dimensionalparameter space. The apparatus of the invention is a registering unitconfigured with the teachings of the method of the invention. Theinvention can also be practiced on a machine readable medium. The methodincludes calculating initial rigid-body transformation parameters for arigid-body transformation on the first image (im₁). The optimumrigid-body transformation parameters are determined by calculating aregistration metric between the second image (im₂) and the rigid-bodytransformation on the first image (im₁). A registered image is generatedfrom the optimum rigid-body transformation parameters.

[0041] In a preferred method for registering 3-d images of a pulmonarynodule from a high-resolution computed tomography (CT) scans, thecalculation of the initial rigid-body transformation parameters ispreferably preceded by masking one of the images by setting pixels to abackground value. Preferably the background value is about −1000. Theregistration metric is generally either minimized or maximized. In onepreferred embodiment, the registration metric is preferably calculatedby:

[0042] transforming the first image (im₁) with the initial rigid-bodytransformation parameters to obtain a transformed first image (im_(1t));

[0043] calculating the registration metric as a correlation (C.) betweenthe transformed first image (im_(1t)) and the second image (im₂); and

[0044] searching for the maximum correlation (C.) in the 6-dimensionalparameter space.

[0045] In another preferred embodiment, the registration metric ispreferably calculated by:

[0046] transforming the first image (im₁) with the initial rigid-bodytransformation parameters to obtain a transformed first image (im_(1t));

[0047] calculating the registration metric as a mean-squared-difference(MSD) between the transformed first image (im_(1t)) and the second image(im₂); and

[0048] searching for the minimum mean-squared-difference (MSD) in the6-dimensional parameter space.

[0049] The transforming of the first image (im₁) to obtain thetransformed first image (im_(1t)) is preferably a mapping of a point vin 3-d space to a point v′ in transformed space defined by:$v^{\prime} = {{R_{x}R_{y}R_{z}v} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}$

[0050] wherein R_(x), R_(y), and R_(z) are rotation matrices defined as:$R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {- {\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}$ $R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{- {\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}$ $R_{z} = {\begin{bmatrix}{\cos \left( r_{z} \right)} & {- {\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}.}$

[0051] The initial rigid-body transformation parameters preferablyinclude six parameters (tx,ty,tz,rx,ry,rz) respectively defined astranslation in x, translation in y, translation in z, rotation about thex-axis, rotation about the y-axis, and rotation about the z-axis.Preferably the initial rotation parameters (rx,ry,rz) are all set tozero, and the initial translation parameters (tx,ty,tz,) are set so thatthe nodule in the first image (im₁) overlaps the nodule in the secondimage (im₂) during the initial calculation of the registration metric.The initial translation parameters (tx,ty,tz,) can be set to adifference between the center of the first image (im₁) and the center ofthe second image (im₂). Preferably the initial translation parameters(tx,ty,tz) are set to a difference between the center of mass of thefirst image (im₁) and the center of mass of the second image (im₂). Thesearching is can be conducted by calculating the correlation (C) ormean-squared-difference (MSD) for every possible set of rigid-bodytransformation parameters. Preferably the searching is conducted byeither a Hill-Climbing search method or by Powell's method.

[0052] The present invention also includes a method and apparatus forremoving extraneous matter from an image having a juxtapleural nodule.The apparatus of the invention is a processing unit configured with theteachings of the method of the invention. The invention can also bepracticed on a machine readable medium. The method includes providing aninitial location P′. A spherical volume that fits inside the image andis centered at the initial location P′ is calculated. A center of mass,COM, of the spherical volume is calculated. An initial direction d′directed towards the extraneous matter is determined in accordance withthe following equation:$d^{\prime} = {\frac{{COM} - P^{\prime}}{{{COM} - P^{\prime}}}.}$

[0053] A current location P_(i) is initialized to be equal to theinitial location P′. A current direction d_(i) is initialized to beequal to the initial direction d′. A maximum ratio γ_(max), step size s,prior mass mass_(i-1), and prior change in mass Δ_(i-1) are initialized.The current location P_(i) is moved by the step size s in the currentdirection d_(i). An equation defining a plane A is determined so thatthe plane A is normal to the current direction d_(i) and the plane Apasses through the current location P_(i). A current mass, mass_(i), iscalculated of the nodule on a side of the plane A opposing that of theextraneous matter. A current change in mass Δ_(i) is calculated bysubtracting the prior mass mass_(i-1) from the current mass mass_(i). Acurrent ratio γ is calculated in accordance with the following equation:$\gamma = {\frac{\Delta_{i}}{\Delta_{i - 1}} - 1.}$

[0054] The prior mass mass_(i-1) is set equal to the current massmass_(i). The prior change in mass Δ_(i-1) is set equal to the currentchange in mass Δ_(i). The current ratio γ is compared to the maximumratio γ_(max). The current direction d_(i) is modified to minimize thecurrent mass mass_(i), and the above steps are repeated starting withmoving the current location P_(i) by the step size s in the currentdirection d_(i) while the current ratio γ is one of less than and equalto the maximum ratio γ_(max). The area of the nodule partitioned by theplane A is output in response to the current ratio γ being greater thanthe maximum ratio γ_(max).

[0055] In the preferred method for removing extraneous matter from animage having a juxtapleural nodule, the extraneous matter can include apleural surface. Preferably the maximum ratio γ_(max) is initialized to0.5 and initializing the step size s is initialized to 1.5. Preferablythe method also includes following the steps after the current ratio γis determined to be greater than the maximum ratio γ_(max):

[0056] defining the current location P_(i) as being visited;

[0057] determining on which side of the plane A the current locationP_(i) is located;

[0058] terminating in reponse to the current location P_(i) not beinglocated on a side of the plane opposing that of the extraneous matter;

[0059] defining the current location P_(i) as being part of a region ofinterest in response to the current location P_(i) being located on theside of the plane opposing that of the extraneous matter; and

[0060] performing these additional steps recursively using a locationcorresponding to at least one of six (6) one-pixel moves from thecurrent location P_(i).

[0061] Preferably the step where the current direction d, is modified tominimize the current mass mass_(i), and the step where the area of thenodule partitioned by the plane A is output in response to the currentratio γ being greater than the maximum ratio γ_(max) include the stepsof:

[0062] calculating recursively the current mass mass_(i) of the noduleon a side of the plane A opposing that of the extraneous matter using atleast one of six (6) directions and a step size s₁ from the currentlocation P_(i); and

[0063] defining the current direction d_(i) equal to the directionyielding the largest decrease in the current mass mass_(i).

[0064] Preferably the initial location P′ is located near a center ofthe nodule.

[0065] For a better understanding of the present invention, reference ismade to the following description to be taken in conjunction with theaccompanying drawings and its scope will be pointed out in the appendedclaims.

BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Preferred embodiments of the invention have been chosen forpurposes of illustration and description and are shown in theaccompanying drawings, wherein:

[0067]FIG. 1 illustrates at the top six (6) images that are selectedslices from a, single computer tomography scan along with thecorresponding extracted lung regions from the selected slices at thebottom;

[0068]FIG. 2 illustrates a histogram of image intensity of the lungregion that includes parenchyma, vessels, and nodules and other solidstructures such as bone and organs;

[0069]FIG. 3 is a table illustrating the densities of structures withinthe lungs;

[0070]FIG. 4 illustrates an air pocket that is not part of the lungs buthas similar characteristics to the modeling of the lungs;

[0071]FIG. 5 is a lung mask generation algorithm;

[0072]FIG. 6 is a flow chart of the lung segmentation algorithm withsample images incorporated therein;

[0073]FIG. 7 illustrates on the left hand side the horizontal streakingartifacts due to beam hardening that can be reduced as shown on theright hand side with median and mean filtering;

[0074]FIG. 8 illustrates an original computer tomography scan on theleft hand side with the image shown after thresholding in the center andthe resulting image after the background has been identified, removedand inverted shown on the right hand side;

[0075]FIG. 9 illustrates on the left hand side a lung after thresholdingwithout median and mean filtering and on the right hand side the imageafter median and mean filtering;

[0076]FIG. 10 illustrates two lungs, gas, and some noise, and that theretention of only the largest components would eliminate the unwantedregions;

[0077]FIG. 11 illustrates at the left hand side a segmented lung afterthresholding but before morphological filtering, and at the right handside the lung after morphological closing;

[0078]FIG. 12 illustrates the lungs being divided medially into threeregions where large vessels are found in region 2 which require a largestructuring element, and that only small vessels are found in regions 1and 3;

[0079]FIG. 13 illustrates a table illustrating the timed execution ofmorphological filtering using both fixed (28 pixel diameter) and varying(28 pixel and 14 pixel diameter) structuring elements;

[0080]FIG. 14 illustrates a matrix containing the histograms of thesegmented lung slices;

[0081]FIG. 15 illustrates an algorithm for finding the location and sizeof a nodule;

[0082]FIG. 16 illustrates a one dimensional model of a) isolated nodulenext to blood vessels, and b) a nodule attached to the pleural wall;

[0083]FIG. 17 illustrates locator template functions;

[0084]FIG. 18 illustrates one dimensional examples of the maximizinglocations of template functions for (a-c) a Gaussian template on anisolated nodule, (d-f) LOG template on a nodule on the wall, and (g-i)LOG template on an isolated nodule;

[0085]FIG. 19 illustrates two dimensional examples of the maximizinglocations for two dimensional LOG templates of various radii on isolatedand pleural nodules where the inner and outer circles represent thepositive and negative regions, respectively, of the LOG template;

[0086]FIG. 20 illustrates two dimensional examples of a sizing templatefunction of various radii which exhibit dramatic changes in thedistribution of values inside the circle when the sizing templatefunction becomes too large;

[0087]FIG. 21 illustrates a detailed algorithm for finding the locationand size of a nodule;

[0088]FIG. 22 illustrates a response of the mean filter function forvarious sized spheres (solid intensity of 1000) versus the radius of thetemplate function;

[0089]FIG. 23 illustrates a response of a Gaussian template function forvarious sized spheres (solid intensity of 1000) versus the radius of thetemplate function;

[0090]FIG. 24 illustrates a graph of the Gaussian and the Laplacian ofthe Gaussian;

[0091]FIG. 25 illustrates a hill climbing search algorithm for thenodule finding algorithm;

[0092]FIG. 26 illustrates locator template functions;

[0093]FIG. 27 illustrates an algorithm for registering a first image anda second image using a rigid-body transformation;

[0094]FIG. 28 illustrates a hill climbing search algorithm fordetermining the best rigid-body transformation parameters for theregistering algorithm;

[0095]FIG. 29 illustrates an algorithm for removing the pleural surfacefrom juxtapleural nodules;

[0096]FIG. 30 illustrates an example of the pleural wall removalalgorithm: The cut nodule is shown as dark gray white the rest of thenodule and the pleural surface are shown in light gray. (A) Given theinitial location and direction, (B) the plane is moved in the direction.(C) The plane eventually intersects the pleural wall causing an increasein the size change of the cut nodule, and (D) the plane is reorientatedto minimize the size of the cut nodule. At the next iteration, (E) theplane intersects the pleural wall, but this time (F) reorientating theplane to minimize the cut nodule still results in a large increase insize change. The algorithm terminates, returning the cut nodule from theprevious iteration, (D);

[0097]FIG. 31 illustrates an algorithm for recursively finding the cutnodule region; and

[0098]FIG. 32 illustrates a hill climbing search algorithm for thealgorithm for removing the pleural surface.

DETAILED DESCRIPTION OF THE INVENTION

[0099] A system in accordance with the present invention may include ascanner, processor, memory, display device, input devices, such as amouse and keyboard, and a bus connecting the various componentstogether. The system may be coupled to a communication medium, such as amodem connected to a phone line, wireless network, or the Internet.

[0100] The present invention is preferably implemented using a generalpurpose digital computer, microprocessor, microcontroller, or digitalsignal processor programmed in accordance with the teachings of thepresent specification, as will be apparent to those skilled in thecomputer art. Appropriate software coding may be readily be prepared byskilled programmers based on the teachings of the present disclosure, aswill be apparent to those skilled in the software art.

[0101] The present invention preferably includes a computer programproduct, which includes a storage medium comprising instructions thatcan be used to direct a computer to perform processes in accordance withthe invention. The storage medium preferably includes, but is notlimited to, any type of disk including floppy disks, optical datacarriers, compact discs (CD), digital video discs (DVD), magneto-opticaldisks, read only memory (ROM), random access memory (RAM), electicallyprogrammable read only memory (EPROM), electrically eraseableprogrammable read only memory (EEPROM), magnetic or optical cards, orany type of media suitable for storing information.

[0102] Stored on any one of the above described storage media, thepresent invention preferably includes programming for controlling boththe hardware of the computer and enabling the computer to interact witha human user. Such programming may include, but is not limited to,software for implementation of device drivers, operating systems, anduser applications. Such storage media preferably further includesprogramming or software instructions to direct the general purposecomputer to perform tasks in accordance with the present invention.

[0103] The programming of the computer preferably includes software fordigitizing and storing images obtained from the image acquisition device(helical computed tomography scanner). Alternatively, it should beunderstood that the present invention may also be implemented to processdigital data derived from images obtained by other means, such as x-raysand magnetic resonance imaging (MRI), positron emission tomography(PET), ultrasound, optical tomography, and electrical impedancetomography.

[0104] The invention may also be implemented by the preparation ofapplication specific integrated circuits (ASIC), field programmable gatearrays (FPGA), or by interconnecting the appropriate component devices,circuits, or modules, as will be apparent to those skilled in the art.

[0105] A. Lung Segmentation

[0106] This section discusses a lung segmentation algorithm in detail.First the model of the lungs will be presented with a discussion of thedifficulties with this task. Then each stage of the algorithm ispresented.

[0107] Referring now to FIG. 1, the object of the algorithm isillustrated. The top 6 images in FIG. 1 are selected slices from asingle CT scan, and the bottom 6 images are the lung regions extractedfrom those slices. The lungs are modeled as a low density regionsurrounded by a high density one. Image filtering is preferably firstperformed to minimize the effects of image noise. A simple lineardiscriminant function based on CT voxel values (photon density) ispreferably used to partition the scan into lung parenchyma and solidstructures. Further filtering can be used to compensate forimperfections in the thresholding.

[0108] The photon density of the structures within the lungs may befound by manually segmenting the lung into different regions andanalyzing the results. Referring to FIG. 3, a table presents the rangeand mean values for the intensity of solid tissue and lung parenchyma.These values were obtained by manually selecting regions within a lungCT scan. Regions of at least 100,000 voxels were selected from four 2.5mm full lung scans. The distribution of intensities is shown in FIG. 2.Note that there is an overlap between the intensity of the lung regionand the intensity of the surrounding solid tissue.

[0109] Referring again to the table in FIG. 3, it is almost possible topartition lung parenchyma and bone using a linear discriminant function.All voxels with a gray level value less than 500 can be consideredparenchyma because the minimum intensity for solid tissue is 545.Unfortunately, vessels, nodules, and partial voxels can not bepartitioned so distinctly. Using the linear discriminant function asdescribed will result in classification errors where the histograms inFIG. 2 overlap, most prominantly between gray levels of 550 to 800. Allvalues of gray level and intensity in the application are in GE unitswhich equates to Hounsfield units HU as follows: HU=GE units −1024, forexample 1000 GE=−24 HU.

[0110] One major challenge is handling the streaking artifacts found in2.5 mm scans with a low radiation dose and a high pitch, such as thosedone in screening studies. When an X-ray beam passes through a regionwith a high attenuation coefficient, the remaining beam contains agreater proportion of high energy. This makes the path-integral ofattenuation a non-linear function of distance (J. Hsieh, “AdaptiveTrimmed Mean Filter for CT Imaging,” SPIE, Vol 2299, 1994; H.Soltanian-Zadeh, J. P Windham and J. Soltanianzadeh, “CT ArtifactCorrection: An Image Processing Approach,” SPIE, Vol 2710, 1996 which isincorporated herein by reference). In the case of lung CT scans, thisartifact is prevalent through the shoulders, resulting in noise in thefirst 25% of the slices. The noise is characterized as both highintensity and low intensity horizontal streaks running through theimage.

[0111] Finally, there are other structures which confound thesegmentation scheme. Air or gas found outside of the lungs will have thesame characteristics as the lung themselves. Specifically they consistof a low density region surrounded by a high density one. It isnecessary to discriminate between the lungs and other stray tissue orair pockets. An example of such a region is shown in FIG. 4, whichoccurs well below the lungs and the diaphragm.

[0112] The segmentation algorithm consists of creating a maskrepresenting the lung volume and applying an AND operation between themask and the original image. The algorithm for creating the mask is setforth in FIG. 5, and FIG. 6 outlines these steps graphically while themotivation behind this routine is described in detail in the followingsections.

[0113] The main object of image filtering is to remove the streakingartifacts on the first 25% of the slices due to x-ray beam hardeningthrough the shoulders. The streaking artifacts are characteristicallyhorizontal and tend to greatly distort the segmentation. Many of thelater steps in this algorithm operate on binary images and are quitesusceptible to noise. A combination of both median and mean filteringreduces the noise significantly. Median filtering preserves the edgesaround the chest cavity while eliminating much of the streaking effect.Since the streaks are horizontal, a small vertical mean filter will blurthe streaking regions, further reducing their effect. This filteringstage is used only to generate a clean image for the mask as theoriginal pixel values are retained in the final segmented lung image.

[0114] Experiments have shown that a median filter size of 4×4 producethe best results. Smaller filters did not reliably eliminate noisecaused by the streaking artifacts. Further experiments have shown that amean filter size of 1×3 produces excellent results. The results offiltering are shown in FIG. 9, but the importance of this step is seenmore clearly in FIG. 7 after the thresholding operation. The filteringis preferably performed on 2 dimensional slices rather than 3D volumes.

[0115] On large data sets, such as full chest CT scans, median filteringis a time consuming operation. To speed up the algorithm, it isimportant to note that the streaking artifacts appears only at thebeginning of the scans. Preferably the median filter for the first 25%of the slices is only computed to reduce the computation time of thisstage by 75%. The mean filter is much smaller (1×3) and thecomputational gains of only running it on part of the image areinsignificant.

[0116] Most of the segmentation is performed by a simple thresholdoperation. As seen in FIG. 2, the distribution of intensities of lungparenchyma and solid tissue is bimodal. A simple threshold of 500 cantherefore be applied to the entire image. For the most part, this willseparate the image into 2 regions: (1) solid tissue and (2) everythingelse. Included in “everything else” are the lungs, trachea, airsurrounding the body, and black background outside of the field of view.By filling in all pixels which touch the image border and inverting theimage, the lungs are clearly segmented from the chest cavity, as shownin FIG. 8.

[0117] As seen in FIG. 2 and the table in FIG. 3, there is some overlapbetween the intensities of structures in the lungs and the intensitiesof solid tissue outside the lungs. Specifically, there is significantoverlap in the range between gray levels of 550 and 800. It is thereforeimpossible to threshold out the ribcage and keep all the vessels, asseen in FIG. 8.

[0118]FIG. 9 illustrates the effect of thresholding both with andwithout image filtering. The noisy image on the left shows the result ofthresholding without filtering and the segmented image on the rightshows the result of thresholding after filtering.

[0119] Up until this stage in the algorithm, any regions with similarcharacteristics to the lungs will be preserved. This includes air andgas found within the body. This also includes random noise in thebackground (i.e., outside the thorax) which was not removed by imagefiltering and thresholding. Selecting only the largest component(s) inthe image, we will eliminate minor air pockets and all stray noise whichis prevalent throughout the images.

[0120]FIG. 10 shows a region which is removed by retaining only thelargest component(s) and discarding all other structures. To perform theselection, preferably all regions with volumes greater than 1% of thetotal number of pixels are chosen in the image, as explained in S. Hu,E. A. Hoffman, and J. M. Reinhardt, “Automatic Lung Segmentation forAccurate Quantitation of Volumetric X-Ray CT Images,” IEEE Transactionson Medical Imaging, Vol 20, No 6, June 2001, which is incorporatedherein by reference.

[0121] As previously mentioned, the thresholding operation will segmentout blood vessels along with the ribcage, resulting in holes andsplotches in the lung region, as seen in FIG. 8. These can be filled inusing a morphological filtering operation. By performing a closing,which consists of a dilation followed by an erosion, the spotty regionsand holes within the lungs can be removed. Some of the vessels and otherstructures are quite large, necessitating the use of a large closingkernel. The results of morphological filtering are shown in FIG. 11.

[0122] Morphological filtering with large structuring elements iscomputationally expensive. Fortunately, large structures that require alarge filter kernel only occur in the middle third of the slices of thelungs. Therefore, the lungs can be divided axially into three distinctregions 1, 2, and 3 as shown in FIG. 12.

[0123] As the slices corresponding to each region show, preferably asmall structuring element is used in regions 1 and 3 where only smallvessels occur, and a larger one is used in region 2, where large vesselsare found. Specifically, a large 2D circular filter is preferably usedin region 2 (typically 28 pixels in diameter) and filter of half thediameter is used in regions 1 and 3. Referring now to FIG. 13, a tablethat shows timing trials using a varying kernel size versus using afixed large kernel for the entire image. In this experiment, a circularfilter of diameter 28 was used as the fixed kernel size. This wascompared to using the same kernel in region two and reducing it to adiameter of 14 in regions 1 and 3.

[0124] Once the binary mask of the lungs is generated, a simple ANDoperation with the original CT scan will yield the segmented lungs.Finally, the image is clipped such that the new bounding box is justlarge enough to contain the lungs.

[0125] B. Lung Volume Measurements

[0126] This section discusses the algorithms used for lung volumemeasurements. First the motivation behind measuring the lung volume willbe explained. Next, a model of the lung will be presented, followed bythe algorithm for calculating the volume.

[0127] Once the lung has been properly segmented, the volume of bothtissue and air in the lungs can be calculated. The objectives of thisresearch are as follows:

[0128] 1. Accurately measure lung volume.

[0129] 2. Determine both the inter-patient and intra-patient variationin lung volume.

[0130] 3. Explore the effects of inspiration on CAD algorithms.

[0131] 4. Compensate for any effects discovered in (3).

[0132] 5. Measure the relationship between cross-sectional area of thelungs and entire lung volume.

[0133] Once statistics of the lung volumes have been compiled the effectof inspiration on nodule volume measurements will be studied. A positivecorrelation between lung volume and nodule size in a large set of repeatscan will determine if such a relationship exists. It is also importantto note that full chest scans are not always performed during a repeatexamination. Typically, a high resolution scan of the suspicious regionis performed. Such scans contain only a few cross-sectional slices ofthe region of interest. Therefore, the relationship between 2Dcross-sectional area and inspiration must be explored.

[0134] The lungs can be modeled as containing air with a gray levelvalue of zero, and tissue with a gray level value of 1000. All intensityvalues in between that range must be a combination of air and tissue,confounded by the partial voxel effect. Based on this assumption, a graylevel of 500 represents a volume comprised of ½ air and ½ tissue.Equation 1 depicts this formally, where V is volume and I is intensity.$\begin{matrix}{V_{{air}/{voxel}} = {\left( {1 - \frac{I_{voxel}}{I_{tissue}}} \right)*\left( V_{voxel} \right)}} & (1)\end{matrix}$

[0135] For example, a gray level of 250 represents 25% tissue and 75%air. If the image resolution is 0.5 mm×0.5 mm×2.5 mm per voxel, thevolume of a voxel is 0.625 mm³ and the volume of air in that voxel is0.156 mm³.

[0136] Referring now to FIG. 14, a matrix is generated from thesegmented lung image in which each column represents a gray level valueand each row represents a slice in the scan to perform this analysis.The entries in the matrix are the number of times that a particular graylevel appears in that slice of the image. This data structure is simplya histogram of the image, separated into slices. With such a structurein place, it is easy to perform and repeat the volume measurements onboth the full lung volume and on an individual slices.

[0137] To calculate the volume in the lungs, the number of voxels in thesegmented image is counted and multiplied by the volume of each voxel.This represents the volume of the entire lung. This can be performedquickly using the histogram because the counting has already beencompleted. Using the model described above, that each voxel is part airand part tissue, the sum of “tissue voxels” multiplied by the volume ofa voxel equals the amount of tissue in the lung. The equations forcalculating the volume of the lungs and the volume of the tissue in thelungs are as follows: $\begin{matrix}{V_{Lung} = {\underset{s \in {Si} \in I}{\sum\sum}{{L\lbrack s\rbrack}\lbrack i\rbrack}*V_{voxel}}} & (2) \\{V_{tissue} = {\underset{s \in {Si} \in I}{\sum\sum}\frac{i}{I_{tissue}}*{{L\lbrack s\rbrack}\lbrack i\rbrack}*V_{voxel}}} & (3)\end{matrix}$

[0138] where V represents volume, I represents voxel intensity, andL[s][i] is the number of pixels of intensity i in slice s of image L.

[0139] Often times with repeat scans, only a high resolution scan of theregion of interest is acquired. As a result, full lung volumes are notalways available to measure the amount of inspiration in repeat scans.As a solution, the area of single axial CT images from two scans may beused to estimate the variation in inspiration between the two. Thehistogram data structure described above can be used to calculate thecross-sectional area of a single slice. Given the volume of the 2 entirelungs, it may be possible to relate the cross-sectional area of each tothe volume of the entire scan.

[0140] The lung segmentation algorithm was run on 196 individual chestscans and 24 pairs of repeat scans from a Cornell University ELCAPdatabase. Out of the 196 scans, 115 contained 2.5 mm slices and 81contained 5 mm slices. The 24 pairs of repeat scans consisted of 19 5 mmscans and 5 2.5 mm scans. Repeat scans of different reconstructions werenot considered. In-plane resolution varied between 0.54 mm/pixel and0.79 mm/pixel.

[0141] Images were acquired on a either a General Electric® LightSpeed™or HiSpeed™ Helical CT scanner. X-ray tube current ranged from 40-200 mAand the tube potential was either 120 or 140 kVp. Testing is preferablyperformed on a dual-processor 700 MHz computer. The algorithm ran inapproximately 6-8 minutes on a full lung scan with 2.5 mm contiguousslices. Without using the varying filter sizes as described above, thealgorithm ran in about 14-15 minutes.

[0142] C. Finding Location and Size of Pulmonary Nodule

[0143] Given an initial seed point inside the nodule, an algorithm findsthe center location of the nodule and the radius of the spherecircumscribing the nodule. The following subsections describe thealgorithm and choices of functions.

[0144] The algorithm is provided a candidate location that is some seedpoint P, which is within the nodule in the image, and a starting radiusr and delivers the optimal location P′ and radius r′ of the nodule. Thealgorithm for finding the location P′ and radius r′ of the nodule isshown in FIG. 15. The value of r is incrementally increased and for eachvalue of r the best location of the nodule center P is computed. A newlocation is calculated by maximizing the correlation (or the response)between the image and a locater template function. The template functionis chosen so that the maximum response is given when the template iscentered over a nodule with radius r as discussed further below. Thus,maximizing the correlation yields the most likely location of a nodulewith the current radius.

[0145] The algorithm determines the radius of the nodule by terminatingitself at the appropriate iteration. Another template function, centeredat P, is applied to the image at the end of each iteration. The sizingtemplate function is designed so that its value decreases dramaticallywhen the radius of the function has become larger than the radius of theactual nodule. Thus, when there is a large decrease in the value of thesecond template function, the radius of the actual nodule has beenexceeded and the algorithm terminates and returns the location andradius from the previous iteration.

[0146] The nodule finding algorithm was developed assuming a two levelmodel where the nodule, blood vessels, and pleural wall are all a singlehigh intensity, while everything else is a single low intensity. Thetemplate functions use parameters of location and radius. The radiusdoes not refer directly to the size of the template function. Rather,the radius is a parameter that describes the size of the nodule forwhich the template will detect; meaning that the template does notnecessarily need to be zero past the extent of the radius

[0147] The locater template functions is designed to locate theapproximate center of mass of a nodule (for nodules attached to thepleural surface a repeatable central location is obtained since thelocation of the actual center of mass is not knowable). The function hasfour parameters: the center in x, y, and z, and the radius. In thealgorithm, the radius is fixed at each iteration, but the center of thetemplate can move. The locater template function was designed, with agiven radius, so that it is maximized when the function is centered overa nodule of the same radius.

[0148] Consider a 1-D example where an isolated nodule is modeled as alarge region of high intensity. FIG. 16A shows an example model wherethe nodule is the center shaded region, while the outer regions consistof noise (composed of high-intensity objects like vessels).

[0149] Consider some simple 1-D template functions shown in FIG. 17. Thesquare template function, FIG. 17A, biases everything equally, causingthe template to be very sensitive to noise in the periphery. Thesensitivity to noise can be reduced by tapering the template function tobe shaped like a triangle, FIG. 17B, or a Gaussian, FIG. 17C. The largevalues near the center will keep the function centered over thehigh-intensity region of the nodule, while the lower values further fromthe center will diminish the effect of any other high-intensity objects.FIGS. 18A through 18C illustrate where the Gaussian template functionwill be maximized for cases where the size of the function is smaller,the same size, and larger than the size of the nodule.

[0150] Now consider a 1-D model of a nodule attached to the pleuralwall, shown in FIG. 18B. The nodule is the lightly shaded centralregion, the pleural wall is the darkly shaded region on the right, andthe lightly shaded region on the left is some noise (blood vessel).Although shaded differently, note that the nodule and the pleural wallare at the same intensity level. With this model, the Gaussian templatefunction will fail to perform correctly. The template does not preventitself from slipping into the wall and centering itself over thecombination of the nodule and the wall.

[0151] Nodules on the pleural wall can be located by choosing a functionthat forces itself away from the wall. This can be achieved by using afunction that has positive values within the radius of the function andnegative values outside the radius, like FIGS. 17D through 17F. Twoexamples of this type of function are the Laplacian of the Gaussian(LOG) and the difference of Gaussians (DOG). When maximizing theresponse of the template function, there are two forces fighting againsteach other. The positive portion of the template function will try tokeep the function centered over the high-intensity regions. Counteringthis, the negative portion of the template will try to push the templateaway from the high-intensity regions, resulting in an equilibrium alongthe edge of the high-intensity region of the nodule. FIGS. 18D through18F show how the LOG-like template function should behave on the 1-Dmodel of a nodule on the wall. If the total weights, the distribution,and the sizes of the positive and negative regions are balancedcorrectly, then the LOG-like template function should behaveaccordingly.

[0152] The LOG-like template function can also be used on isolatednodules. The positive region of the template function will stay near thecenter of the nodule, while the negative regions will force the functionaway from any noise in the periphery of the image. FIGS. 18G through 18Ishow how the LOG-like function will behave when applied to the 1D modelof an isolated nodule. When the size of the template is smaller than thenodule, the template will hug the edge of the nodule and when thepositive region of the template is larger than the nodule, it will becentered over the nodule.

[0153] The above illustrates the preferred behavior of the locatortemplate function for the 1D model of an isolated nodule and a nodule onthe wall. The model is now extended into two and three dimensions. Forthe 2D case, the 1D LOG-like function is extended into two dimensions.The template now consists of two circular regions; the smaller regioncontains positive values and the region outside the smaller regioncontains negative values. FIG. 19 shows a few examples of how the 2Dtemplate function is expected to behave for isolated and pleuralnodules. For the 3D case, the locater template function will havespherical regions with an analogous distribution as the 2D and 1D cases.

[0154] One thing to note about the locater template function is that itdoes not need to have a constant total weight for different radii. Thealgorithm does not maximize the locater template function throughout allthe iterations and all the radii; instead, it is maximizing the templatefunction for a fixed radius in each iteration. Thus, the template is notrequired to have a constant weight for different radii. It shouldsuffice to only maintain the proper scaling between the sizes and thedistribution of the template regions.

[0155] The sizing template function is used to indicate when the radiusof the function is equal to the size of the sphere circumscribing thereal nodule. The sizing template function uses the radius during thecurrent iteration and the location of the candidate nodule, determinedby the locater template function. After each iteration in the algorithm,a value is calculated using the sizing template function. Two preferredmethods for determining when the algorithm has found the location andsize of the nodule in the image are discussed below.

[0156] The first method is to look for a large change in the response ofthe sizing template function. Using the LOG-like locater templatefunction, the positive region, the inner circle, is expected to beoverlayed on the high-intensity values of the nodule. As the radius ofthe template function increases, the positive region will remain withinthe nodule. This is illustrated in FIG. 20. Eventually the positiveregion will become larger than the nodule and there will be a sharpdecrease in the response of the sizing template function betweeniterations of the algorithm. To find this stopping criterion, a sizingtemplate function is chosen so to be a mean filter shaped like a sphere.When the sphere becomes larger than the nodule, the distribution ofvalues inside the sphere changes, causing the response of the templatefunction to change dramatically.

[0157]FIG. 22 shows the response of the mean filter function to varioussizes of ideal nodules (spheres of intensity 1000). For the idealnodule, the response of the filter function is constant until the radiusof the filter function exceeds the radius of the ideal nodule. When thishappens, the response of the filter quickly decreases.

[0158] The second method of finding the size of a nodule is to designthe sizing template function so that it has a specific response when thetemplate function is the same size as the real nodule. For example, if a3D Gaussian function, where the standard deviation is equal to theradius, has the response shown in FIG. 23. When the response of theGaussian template is 200, the radius of the Gaussian is the same as theradius of the nodule, thus providing a stopping criterion for ouralgorithm.

[0159] One very important property of the sizing template function isthat it must have a constant total weight regardless of the size of thetemplate. This property will insure that the response of the templatewill be consistent between different sizes of nodules.

[0160] The implementation of the algorithm will now be discussed. Thiswill include a more detailed description of the algorithm, amathematical definition of some template functions, calculation of theinitial seed point, introduction of two search methods for themaximization of the locator template function, quick functionevaluations, and the termination criterion.

[0161] A detailed description of the algorithm is shown in FIG. 21. Thealgorithm is modified so that it terminates when the response of thesizing template function is within c of the target value. This featureis added to the algorithm so that the radius does not go too far pastthe actual radius of the nodule.

[0162] The algorithm takes in a 3d CT image of a nodule. First, theimage is clipped at intensities over 1000. This is done to minimize anyaffect that the very-high intensitied bone will have on the algorithm.

[0163] The next step is to define the locator and sizing templatefunctions. The template functions are preferably either the Gaussian,the Laplacian of the Gaussian, or the difference of Gaussians. Thetemplate function has four parameters, (sx, sy, sz, σ), which correspondto the x-location, y-location, z-location, and radius. The initiallocation is either calculated from the image or specified by the user,and depending on the type of sizing template function, the targetresponse value is set.

[0164] Starting with the initial conditions, a radius of 1, and a radiusstep-size of 1, the algorithm finds the best locator template responseby using a search method to move the center of the template function.Preferably the search method is a hill-climbing method as discussedbelow. Once the maximum value is found for the given radius, the radiusis increased by the radius step-size and the search is repeated usingthe current location as the starting point.

[0165] At each iteration the response of the sizing template function, γis calculated. The algorithm iterates until the sizing response is lessthan the target value. If γ is within ε of the target, then thealgorithm is finished. Otherwise, the algorithm backtracks by settingthe location and radius parameters to their values in the previousiteration, dividing the radius step by two, and repeating the searchprocedure. Finally, after the target value is reached within ε, or theradius step-size is smaller than some α, the search process terminatesand the current parameters are the location and radius of the nodule.

[0166] The template function is a function of four parameters: thecenter (location) of the template, (sx,sy,sz), and the radius of thetemplate, σ. The response of the template function is calculated bytaking the correlation between the image im and the template M as inEquation (4). $\begin{matrix}{{\Gamma \left( {{sx},{sy},{sz},\sigma} \right)} = {\frac{\sum{\sum\sum}}{x,y,z}{{im}\left( {x,y,z} \right)}*{M_{{sx},{sy},\sigma}\left( {x,y,z} \right)}}} & (4)\end{matrix}$

[0167] The locater template function is designed to produce the largestresponse when the kernal is centered over the target object. Thus, bymaximizing this function at each iteration of the nodule findingalgorithm, the best candidate for a nodule of radius r is found. Thereare many types of functions to choose from, but the preferred functionsare the Gaussian, the Laplacian of the Gaussian, and the Difference ofGaussians.

[0168] The 3D gaussian is a strictly positive symmetric function definedby Equation (5). The Gaussian function can be used to find nodules whenthe pleural surface is not present in the image. The metric functionwill have its highest value when it is centered over the nodule. A graphof the Gaussian is shown in FIG. 24. $\begin{matrix}{{G_{{sx},{sy},{sz},\sigma}\left( {x,y,z,} \right)} = {\frac{1}{{\sigma^{3\quad}\left( {2\pi} \right)}^{\frac{3}{2}}}^{\frac{- {({{({x - {sx}})}^{2} + {({y - {sy}})}^{2} + {({z - {sz}})}^{2}}}}{2\sigma^{2}}}}} & (5)\end{matrix}$

[0169] The Laplacian of the Gaussian (LOG) is defined in Equation 7.This function has a positive weight close to the center and a negativeweight further out, as seen in FIG. 24. This function is useful inlocating nodules on or near the pleural surface. The positive interiorwill latch onto the nodule, while the negative exterior will keep thekernal from moving towards the wall. In the algorithm, the LOG has thedisadvantage that the initial location must be within the nodule.Otherwise, the negative weight of the exterior will cause the LOG togrow away from the nodule. $\begin{matrix}{{LOG}_{{sx},{sy},{sz},{\sigma {({x,y,z})}}} = {{- \left\lbrack {\frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}} + \frac{\partial^{2}}{\partial z^{2}}} \right\rbrack}{G_{{sx},{sy},{sz},\sigma}\left( {x,y,z,} \right)}}} & (6) \\{= {\frac{1}{\sigma^{2}}\left( {3 - \frac{x^{2}}{\sigma^{2}} - \frac{y^{2}}{\sigma^{2}} - \frac{z^{2}}{\sigma^{2}}} \right)*{G_{{sx},{sy},{sz},\sigma}\left( {x,y,z} \right)}}} & (7)\end{matrix}$

[0170] The Difference of Gaussians, given in Equation (8), is similar tothe Laplacian of the Gaussian in that they both have positive andnegative regions. However, changing the σ₁, σ₂, m₁, and m₂ values willlead to different sizes and weights of the two regions, allowing for amore customizable template than the LOG.

DOG _(sx,sy,sz,σ)(x,y,z)=m ₁ *G _(sx,sy,sz,σ1)(x,y,z)−m ₂ *G_(sx,sy,sz,σ2)(x,y,z)  (8)

[0171] The initial seed point should be somewhere within the nodule inthe image. This may be provided by user input, or it can be estimatedfrom the image. A simple way to calculate the initial location would beto use the center of the image:

[0172] center_(x)=im.xlo+(im.xhi−im.xlo+1)/2

[0173] center_(y)=im.ylo+(im.yhi−im.ylo+1)/2

[0174] center_(z)=im.zlo+(im.zhi−im.zlo+1)/2

[0175] When using the Gaussian for the locator template function thismay be fine, but when using the LOG or DOG functions the initialstarting location must be within the nodule. For these functions, thecenter of mass (COM) of the image is preferred for setting the initialconditions. The COM is calculated by first thresholding the image athalf of the largest pixel value in the image and then calculating thecenter of mass using the standard equations:${COM}_{x} = \frac{\sum{\sum{\sum_{i,j,k}{{{im}\left( {i,j,k} \right)}*i}}}}{\sum{\sum{\sum_{i,j,k}{{im}\left( {i,j,k} \right)}}}}$${COM}_{y} = \frac{\sum{\sum{\sum_{i,j,k}{{{im}\left( {i,j,k} \right)}*j}}}}{\sum{\sum{\sum_{i,j,k}{{im}\left( {i,j,k} \right)}}}}$${COM}_{z} = {\frac{\sum{\sum{\sum_{i,j,k}{{{im}\left( {i,j,k} \right)}*k}}}}{\sum{\sum{\sum_{i,j,k}{{im}\left( {i,j,k} \right)}}}}.}$

[0176] The method used to search for the optimal nodule location has amajor impact on the running time of the algorithm. Any kind of fancysearch procedure could be used, but this is not necessary for thisapplication. From iteration to iteration the new optimal location willnot deviate much from the previous optimal location. Since the previousoptimal location is used as the starting point for each search, thesearch method will not have to look very far for the maximum in thesearch-space.

[0177] Fancy search methods, such as Powell's method described in W.Press, Numerical Recipes in C, 2nd Edition, Cambridge University Press,1992, which is incorporated herein by reference, are usually lessefficient when moving small distances because these algorithms aredesigned to search a very large space. In other words, they areoptimized to tradeoff between the efficiency in moving giant steps andthe efficiency of moving in small steps.

[0178] In the end, hill-climbing, a greedy search method, is preferredfor the nodule finding algorithm because it is more efficient when theoptimal location doesn't move very far. The Hill-Climbing is shown inFIG. 25. Using an the initial location and a stepsize, the hill-climbingalgorithm evaluates the locator template function for each of thepossible 6 moves in parameter-space and then moves in the direction ofthe largest decrease. The algorithm terminates when there is no movethat decreases the metric.

[0179] The search procedure requires many evaluations of the correlationbetween locator template function and the image. These calculations canbe costly because the template function usually involve exponentials.Fortunately, in any given iteration of the search procedure, the radiusof the template function is held constant. It is possible to savecalculations by storing the template function values in a temporaryimage (the classic engineering dichotomy of time versus space).

[0180] The values of the locater template function centered at theorigin are stored in a temporary image at the beginning of each searchprocedure. For the Gaussian and the LOG, values are only calculated outto 3 or 5 times the radius σ because the functions converge to zero.Also, the function is only evaluated at positive location points becausethese functions are radially symmetric,

[0181] Now when the response of the template function is required, thetemporary image is offset by the location of the template function, andthe correlation between the image and the temporary image is calculated.

[0182] The target response determines when the algorithm shouldterminate. The target value is chosen such that the algorithm willterminate when the location and radius parameters best circumscribes thenodule. The target value can be found either through an ideal model of anodule or through experimentation. FIG. 23 shows the response of aGaussian function to several ideal nodules (spheres) of different radii.In the ideal case, the target response of the Gaussian sizing functionshould be 200. The best target response has been determined to be 300 byexperimentation with a small subset of real nodules from the database.Using a small number (four) of nodules, the response of the sizingfunction was tracked and the value that causes all the nodules to becircumscribed was selected.

[0183] Based upon the experimentation, the locator template functionsshown in FIG. 26 will be stable. The large negative values near theborder between the positive and negative regions help the template findthe edges of the nodule. In addition, the difference of Gaussians (DOG)can be used to correctly balance the weights and sizes of the positiveand negative regions in the template so that the template will not slipinto the pleural wall

[0184] The nodule finding algorithm is preferably implemented in the Cprogramming language for a VisionX software package on a FreeBSD system.VisionX software provides computer tools and programs for the analysisand visualization of image data. It is suitable for a wide range ofimage analysis applications and is designed to address the processingneeds of multidimensional image sets that arise both from temporal imagesequences and from image modalities that involve three-dimensional datacollection.

[0185] VisionX software has been used in a wide range of researchapplications including multispectral image analysis, three-dimensionalobject recognition, multiframe image analysis, target tracking, neuralnetworks, biological cell analysis, and three-dimensional biomedicalimage analysis. Important features of the VisionX software include theability to handle multidimensional image sets, a wide range of availableprocessing functions, and a flexible tagged data format that facilitiesthe automatic recording of the history of a file.

[0186] FreeBSD is an advanced operating system for x86 compatible, DECAlpha, and PC-98 architectures. It is derived from BSD UNIX, which is aversion of UNIX developed at the University of California, Berkeley.FreeBSD offers advanced networking, performance, security, andcompatibility features. FreeBSD can be installed from a variety of mediaincluding CD-ROM, DVD-ROM, floppy disk, magnetic tape, an MS-DOSpartition, or if there is a network connection, it can be installeddirectly over anonymous FTP or NFS.

[0187] The nodule finding algorithm preferably inputs a floating point,byte, or short image and outputs the location and size of the nodule.Many options are available, including translation and radius searchlimits, a choice between search strategies (Hill-Climbing or Powell'smethod), initial starting location, type of template function, andtermination criteria. The program also preferably outputs several typesof images for debugging purposes. The program was expanded to allow thealgorithm to run on two-dimensional images.

[0188] D. Pulmonary Nodule Registration

[0189] The algorithm registers pulmonary nodules in 3-d images ofhigh-resolution focused CT-scans. For the most accurate characterizationit is important to have images from high-resolution CT scans becausethese images will better localize the boundaries of the nodule. Arigid-body transformation model is assumed, meaning that in general thestructures in the images are confined to simple translation androtation. In order to simplify the transformation model, an isotropicimage-space is also assumed. Using these assumptions, an algorithm wasdeveloped to register pulmonary nodules from two different time periodsby defining a metric between the two images and performing a minimizingsearch on the metric.

[0190] The algorithm requires two input images, and produces an outputthat is the first image registered to the second image. The rigid-bodyregistration algorithm is shown in FIG. 27. The two input images, im₁and im₂, are converted to floating-point pixel format. Preferably theimages are masked to ignore irrelevant pixel data (e.g. bones). Next,the initial conditions for the rigid-body transformation is determined.The metric between two images is defined as a number that represents howclosely two images are related to each other. By selecting theappropriate rigid-body parameters, the metric between the second imageand the transformation of first image can be minimized or maximizeddepending upon the registration metric, resulting in the bestregistration for the given metric.

[0191] The first step in the registration algorithm is to perform somepre-processing on the input images. If necessary, the input image isconverted from byte or short pixel-formats to the floating pointpixel-format.

[0192] If mask images are provided, the input images are masked bysetting the appropriate pixels to the background value (−1000). Maskingan image tells the registration algorithm to ignore certain areas ofeither image when attempting to perform the registration. This can beuseful if a large structure in an image is causing misregistration ofthe object of interest. By using a mask, the algorithm can be told toignore the larger structure, thus allowing the object of interest to beregistered to itself.

[0193] The registration metric is defined as a function of two images,im₁ and im₂ and six rigid-body transformation parameters,(t_(x),t_(y),t_(z),r_(x),r_(y),r_(z)). To calculate the metric, thefirst image is transformed using the rigid-body parameters, resulting inim_(1t). One definition of the metric value is the correlation betweenthe transformed image im_(1t) and the second image im₂, as given in thefollowing equation. $\begin{matrix}{C = {\frac{1}{N}\begin{matrix}{\sum{\sum\limits_{i,j,k}\quad {\sum{i\quad {m_{1t}\left( {i,j,k} \right)}*i\quad {m_{2}\left( {i,j,k} \right)}}}}} \\{{i\quad {m_{1t}\left( {i,j,k} \right)}} \neq {background}} \\{{i\quad {m_{2t}\left( {i,j,k} \right)}} \neq {background}}\end{matrix}}} & (9)\end{matrix}$

[0194] where N is the number of pixels that are not background pixels ineither im_(1t) or im₂.

[0195] The correlation metric is not absolute in the sense that there isno absolute best correlation that indicates that the two images areregistered perfectly. The correlation is dependant on the distributionof the pixels; images with larger pixel values will produce a largercorrelation. An absolute metric allows for a better sense of how welltwo images are registered to each other regardless of pixeldistributions.

[0196] Instead of correlation, a mean-squared-difference (MSD) metric,defined in the following equation, can be used: $\begin{matrix}{{MSD} = {\frac{1}{N}\begin{matrix}{\sum{\sum\limits_{i,j,k}\quad {\sum\left( {{i\quad {m_{1t}\left( {i,j,k} \right)}} - {i\quad {m_{2}\left( {i,j,k} \right)}}} \right)^{2}}}} \\{{i\quad {m_{1t}\left( {i,j,k} \right)}} \neq {background}} \\{{i\quad {m_{2t}\left( {i,j,k} \right)}} \neq {background}}\end{matrix}}} & (10)\end{matrix}$

[0197] where N is the number of pixels that are not background pixels ineither im_(1t) or im₂. With the MSD metric, perfectly registered imageswill produce a metric of zero.

[0198] The rigid-body transformation uses six parameters:(tx,ty,tz,rx,ry,rz); translation in x, translation in y, translation inz, rotation about the x-axis, rotation about the y-axis, and rotationabout the z-axis. The transformation is a mapping of a point v in 3-dspace to a point v′ in transformed space defined by the followingequation: $\begin{matrix}{\nu^{\prime} = {{R_{x}R_{y}R_{z}\nu} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}} & (11)\end{matrix}$

[0199] where R_(x), R_(y), and R_(z) are the rotation matrices definedas: $\begin{matrix}{R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {- {\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}} & (12) \\{R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{- {\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}} & (13) \\{R_{z} = \begin{bmatrix}{\cos \left( r_{z} \right)} & {- {\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}_{.}} & (14)\end{matrix}$

[0200] The first image is transformed so that it has the same boundingbox as the second image. The transformation uses either linearinterpolation or nearest-neighbor interpolation, and pixels that aretransformed from outside the bounding box of the first image are set tothe background value.

[0201] The initial rigid-body transformation parameters must bedetermined before a minimization search can be performed. The initialrotation parameters are all set to zero because only a small amount ofrotation is expected between the two images. The initial translationparameters should be set so that the two nodules will overlap at thefirst iteration of the search procedure. Assuming that the nodules arein the center of the images, the initial translation parameters can beset to the difference between the centers of the two images, where thecenter of an image is defined as:

[0202] center_(x)=im.xlo+(im.xhi−im.xlo+1)/2

[0203] center_(y)=im.ylo+(im.yhi−im.ylo+1)/2

[0204] center_(z)=im.zlo+(im.zhi−im.zlo+1)/2

[0205] This assumption is not always true and a better method ofdetermining the initial translation parameters is to use the centers ofmass of the two nodules. Both images are thresholded at half of thelargest pixel value in the image and the center of mass is calculatedusing the standard equations: $\begin{matrix}{{COM}_{x} = \frac{\sum{\sum{\sum_{i,j,k}{i\quad {m_{\quad}\left( {i,j,k} \right)}*i}}}}{\sum{\sum\quad {\sum_{i,j,k}{i\quad {m_{\quad}\left( {i,j,k} \right)}}}}}} \\{{COM}_{y} = \frac{\sum{\sum{\sum_{i,j,k}{i\quad {m_{\quad}\left( {i,j,k} \right)}*j}}}}{\sum{\sum{\sum_{i,j,k}{i\quad {m_{\quad}\left( {i,j,k} \right)}}}}}} \\{{COM}_{z} = \frac{\sum{\sum{\sum_{i,j,k}{i\quad {m_{\quad}\left( {i,j,k} \right)}*k}}}}{\sum{\sum\quad {\sum_{i,j,k}{i\quad {m_{\quad}\left( {i,j,k} \right)}}}}}}\end{matrix}_{.}$

[0206] The initial translation parameter is equal to the differencebetween the centers of mass of the two images.

[0207] Using the registration metric, a function of six parameters, andthe initial conditions, the registration problem is now converted into aminimization/maximization search problem. Specifically, given someinitial condition, a minimum or maximum of the metric function is foundby searching the 6-dimensional parameter space. The followingsubsections describe three different search methods: exhaustive search,hill-climbing, and Powell's method.

[0208] The exhaustive search is a brute force method that finds theminimum (or maximum) by calculating the metric for every possible set ofparameters. Of course the number of possible parameters in a6-dimensional search space is infinite. In order to reduce the searchspace, a limit of translation and rotation deviation from the initialcondition, along with rotation and translation stepsizes, is imposed.Even with this search space reduction, the exhaustive search is stillimpractical.

[0209] Consider a case where translation is limited to 5 pixels in eachdirection and the rotation by 5 degrees in each direction. Withstepsizes of 1 pixel of translation and 1 degree of rotation, theexhaustive search will perform 11⁶ or 1,771,561 metric evaluations.Suppose each metric evaluation takes 0.1 seconds, then the entirealgorithm will take around 170,000 seconds, or 2,833 minutes, or justunder two days. Even with such a seemingly small search space andappropriate step-sizes, the exhaustive search will take days to finish!

[0210] Hill-Climbing is a greedy search method. Starting at the initialcondition and a set of translation and rotation stepsizes, the algorithmevaluates the metric for each of the possible 12 moves inparameter-space and then moves in the direction of the largest decrease(or increase). The algorithm terminates when there is no move thatimproves the metric.

[0211] By using a small step-size it is possible to calculate the bestparameters up to some precision. Unfortunately, as the step-size isdecreased, the running-time of the algorithm increases because thealgorithm must take smaller steps, and thus perform more metricevaluations to get to the minimum. One solution to this problem is torepeat the algorithm several times while using a decreasing stepsize. Adetailed description of the complete hill-climbing algorithm is shown inFIG. 28.

[0212] Powell's method is a multi-dimensional direction-set searchalgorithm. Starting with an initial set of 6 directions and an initialcondition, in each iteration the algorithm minimizes the metric bymoving in each of the six directions. For a given direction, any lineminimization technique could be used, however Brent's Method is chosenbecause it is a parabolic minimization technique that does notexplicitly use derivatives. At the end of each iteration (one passthrough each of the six directions), the oldest direction is replaced bythe total direction moved during the current iteration. By doing this,the algorithm adapts itself to move along the most minimizing path. Amore in depth description of Powell's method can be found in Chapter10.5 in Numerical Recipes in C, W. Press, 2^(nd) Ed., CambridgeUniversity Press, 1992, which is incorporated herein by reference.

[0213] Powell's method terminates when either all the parameters in oneiteration change within some epsilon, or the metric after an iterationdoes not change by more than a tolerance value.

[0214] The registration algorithm works well on images with isolatednodules, but has trouble with images where the pleural surface ispresent because the pleural surface becomes the largest structure in theimage. The largest structure has the largest effect on the metric value,causing the minimization technique to try to register the pleuralsurface to itself. Because the pleural wall can move and change shapewith inspiration or position, there is no guarantee that the nodulesnear or along the pleural wall be registered correctly if the pleuralwall is registered. A first method for addressing this problem is to useimage masking or pixel masking. A mask is created that tells thealgorithm to ignore the pleural surface and its contents will force thealgorithm to register only the nodules together. A second methodinvolves applying a function on the pixels so that all the pleuralfeatures become one intensity. This causes the registration algorithm toignore structures like the ribs. A third method involves reducing thesearch space. A nodule-localization algorithm is used to find thecenters of the nodules in both images. Taking the difference of thecenters results in the translation parameters that overlap the nodules.Finally, the orientation of the nodules is determined by using theregistration algorithm on only the rotation parameters.

[0215] The rigid-body registration algorithm is preferably implementedin the C programming language using the VisionX software library for theFreeBSD environment. The program, v3regrb, preferably inputs two inputimages of byte, short, or floating point pixel-types and a variety ofoptions. The program preferably outputs the first input image registeredto the second input image and a difference image associated with the twoimages. Both images preferably have the same bounding box as the secondinput image. The rigid-body transformation parameters for theregistration are preferably stored with the metric and timinginformation with a history of the output images.

[0216] In order to increase the throughput of the registration process,the three-dimensional rigid-body transformation is preferably writteninto the program. The code for the transformation is preferably usedfrom the existing stand-alone rigid-body transformation program(v3regrb). The rigid-body function preferably supports both linear andnearest-neighbor interpolation. The program may be modified to performthree-dimensional translation-only registration, two-dimensionalrigid-body registration, and two-dimensional translation-onlyregistration.

[0217] E. Removal of the Pleural Surface from Juxtapleural Nodules inThresholded High-Resolution CT Images

[0218] The algorithm takes a binary (thresholded) image as its input.The algorithm also requires the location of a point near the center ofthe nodule. If the image is some region-of-interest that has beengenerated by a radiologist, then it is assumed that the nodule is in thecenter of the image. In this case the initial starting point will be thecenter of the image. The pleural-surface removal algorithm is shown inFIG. 29.

[0219] The pleural surface removal algorithm works by iteratively movinga plane towards the pleural-surface. The algorithm starts with aninitial point P′ inside the nodule and a direction towards the wall, d′.The direction is calculated by taking the difference between the centerof mass of a spherical region centered on P′ and the starting locationP′.

[0220] With each iteration, a new location P_(i) is calculated bystepping in direction d from the previous location P_(i-1). Plane A,normal to direction d and passing through point P, separates the nodulefrom the pleural surface. This cut nodule is the connected componentregion that contains point P and that is behind the plane A. As theplane moves towards the pleural wall, the size of the cut noduleincreases. FIGS. 30A and 30B show the cut nodule for two iterations.

[0221] The algorithm keeps track of the difference Δ in the size of thecut nodule between iterations. When the plane intersects the pleuralsurface, as seen in FIG. 30C, the difference Δ increases dramatically.This condition is manifested by the change in the difference increasingby more than γ_(max). When this occurs, the size of the cut nodule isminimized by reorientating the plane by changing the direction d whilekeeping the point P fixed, as in FIG. 30D. Additionally, the sizedifference and change in difference are recalculated. If the change indifference is less than γ_(max), then the algorithm keeps iterating.Otherwise, the algorithm terminates and returns the cut nodule formed byusing the plane in the previous iteration.

[0222]FIG. 30 shows a two-dimensional example of the algorithm. The cutnodule is shown as dark gray while the rest of the nodule and thepleural surface are light gray. The initial starting point and theinitial direction are shown in FIG. 30A. FIG. 30B shows an iterationwhere the plane does not intersect the pleural wall. FIG. 30C shows aniteration where the plane intersects the pleural surface, causing thechange in the nodule size to increase. In FIG. 30D the plane isreorientated to minimize the cut nodule, forming a new direction. Afteranother iteration, FIG. 30E shows the plane intersects the pleural wallagain, but reorientating the plane still leads to a large increase innodule size, as seen in FIG. 30F. Thus, the algorithm terminates andreturns the cut nodule from the previous iteration, FIG. 30D.

[0223] In this section, a specific implementation of the algorithm isdiscussed. A seeded region growing algorithm is used to find thecut-nodule region, and a hill-climbing search is used to perform theminimization when reorientating the plane. This section first discussesthe equations used for the plane, followed by a description of theregion growing and hill-climbing algorithms.

[0224] The representation of the plane has four parameters and isdefined as:

ax+by+cz+d=0  (15)

[0225] Given a direction normal to the plane, v, and a point on theplane, P, the parameters of the plane are calculated as:

a=v_(x)  (16)

b=v_(y)  (17)

c=v_(z)  (18)

d=−(aP _(x) +bP _(y) +cP _(z))  (19)

[0226] Finally, a point p is behind the plane if:

ap _(x) +bp _(y) +cp _(z) +d<0  (20)

[0227] A recursive region growing is used to determine the cut-noduleregion. The starting point of the region-growing algorithm is theinitial starting location P′, and the plane A is calculated from P_(i)and d_(i). A description of the region growing algorithm is shown inFIG. 31. For the given pixel p, the algorithm marks it as visited andthen determines if it is behind the plane A. If it is not, then thealgorithm exits. Otherwise, the pixel p is marked as part of thecut-nodule region, and the algorithm is recursively called for each ofthe six possible 1-pixel moves from p if the new point has not beenvisited yet and it is also a foreground pixel.

[0228] A simple greedy search algorithm is used to reorientate the planeso that it minimizes the cut-nodule size. A new plane is calculated bychanging the direction normal to the plane. The hill-climbing algorithmis shown in FIG. 32. The hill-climbing algorithm takes the initialdirection and a small stepsize. Next, the size of the cut-nodule iscalculated for the planes formed by moving the direction in the sixpossible coordinate directions. If a smaller size is calculated, thenthe direction is moved in the direction of the largest decrease. The newdirection is normalized to one, and the algorithm is repeated. If thereis no decrease in the cut-nodule size, then the algorithm terminates.

[0229] The algorithm is preferably implemented in the C programminglanguage for the VisionX software package on a FreeBSD UNIX system. Theprogram preferably operates on binary images that have floating point,byte, or short pixel types. The program preferably outputs the segmentednodule and two debugging images. The initial conditions, initialdirection, and termination criteria are preferably specified by options.

[0230] Thus, while there have been described what are presently believedto be the preferred embodiments of the invention, those skilled in theart will realize that changes and modifications may be made theretowithout departing from the spirit of the invention, and is intended toclaim all such changes and modifications as fall within the true scopeof the invention.

What is claimed is:
 1. A method for generating a lung mask forsegmenting a lung image from voxel data containing the lung image,noise, solid components and surrounding background information, themethod comprising: (a) applying a median filter and a mean filter to thevoxel data to reduce noise; (b) thresholding the noise reduced voxeldata to identify solid components and other structures; (c) identifyingthe surrounding background in the noise reduced voxel data; (d) deletingthe surrounding background from the noise reduced voxel data; (e)labeling connected components of the voxel data; (f) determining thelargest connected components to select a lung region from the voxeldata, the lung region having a geometric form; and (g) morphologicalfiltering the voxel data to refine the geometric form of the lungregion.
 2. The method as defined by claim 1, wherein the voxel data isassociated with a plurality of slices through a patient's lung, theslices beginning at about the patient's shoulders; and step (a) isperformed for the first 25 percent of the plurality of slices onlywhereby the computation time is substantially reduced.
 3. The method asdefined by claim 1, wherein the median filter has a size of 4×4.
 4. Themethod as defined by claim 1, wherein the mean filter has a size of 1×3.5. The method as defined by claim 1, wherein the voxel data in step (b)is thresholded at a gray level of about
 500. 6. The method as defined byclaim 1, wherein the largest connected components are associated withmore than about 1 percent of the voxel data.
 7. The method as defined byclaim 1, wherein the voxel data is associated with a plurality ofslices, the slices being divided into a first end region, a middleregion, and a second end region; and the morphological filtering beingperformed with 2D circular filter having a first diameter in the firstend region and the second end region; and the morphological filteringbeing performed with 2D circular filter having a second diameter whichis about twice the first diameter in the middle region.
 8. A method formeasuring lung volume from a segmented lung image obtained from a scanwhich includes a plurality of slices of voxel data, the voxel datahaving a gray level value and a volume associated therewith, the methodcomprising: generating a matrix of entries from the segmented lungimage, the matrix having a plurality of columns and a plurality of rows,each of the plurality of columns representing the gray level value andeach of the plurality of rows representing one of the plurality ofslices in the scan, the entries corresponding to the number of timesthat the gray level occurs in the corresponding slice; determining anumber of voxels in the segmented lung image from the matrix entries;and multiplying the number of voxels in the segmented lung image by thevolume of each voxel.
 9. A method for measuring volume of tissue in asegmented lung image obtained from a scan including a plurality ofslices of voxel data, the voxel data having a gray level value and avolume, the method comprising: generating a matrix of entries from thesegmented lung image, the matrix having a plurality of columns and aplurality of rows, each of the plurality of columns representing thegray level value and each of the plurality of rows representing one ofthe plurality of slices in the scan, the entries corresponding to thenumber of times that the gray level occurs in the corresponding slice;determining a sum of tissue of voxels in the segmented lung image fromthe matrix entries; and multiplying the sum of tissue of voxels in thesegmented lung image by the volume of each voxel.
 10. The method asdefined by claim 9, wherein the sum of tissue of voxels is calculated bysumming the product of each matrix entry multiplied by the correspondinggray level value divided by a gray level value assigned for tissue. 11.A lung mask generating apparatus for segmenting a lung image from voxeldata containing the lung image, noise, solid components and surroundingbackground information, the lung mask generating apparatus comprising: amasking unit configured to: (a) apply a median filter and a mean filterto the voxel data to reduce noise; (b) threshold the noise reduced voxeldata to identify solid components and other structures; (c) identify thesurrounding background in the noise reduced voxel data; (d) delete thesurrounding background from the noise reduced voxel data; (e) labelconnected components of the voxel data; (f) determine the largestconnected components to select a lung region from the voxel data, thelung region having a geometric form; and (g) morphologically filter thevoxel data to refine the geometric form of the lung region.
 12. A lungmask generating apparatus as defined by claim 11, wherein the voxel datais associated with a plurality of slices through a patient's lung, theslices beginning at about the patient's shoulders; and step (a) isperformed for the first 25 percent of the plurality of slices onlywhereby the computation time is substantially reduced.
 13. A lung maskgenerating apparatus as defined by claim 11, wherein the median filterhas a size of 4×4.
 14. A lung mask generating apparatus as defined byclaim 11, wherein the mean filter has a size of 1×3.
 15. A lung maskgenerating apparatus as defined by claim 11, wherein the voxel data instep (b) is thresholded at a gray level of about
 500. 16. A lung maskgenerating apparatus as defined by claim 11, wherein the largestconnected components are associated with more than about 1 percent ofthe voxel data.
 17. A lung mask generating apparatus as defined by claim11, wherein the voxel data is associated with a plurality of slices, theslices being divided into a first end region, a middle region, and asecond end region; and the morphological filtering being performed with2D circular filter having a first diameter in the first end region andthe second end region; and the morphological filtering being performedwith 2D circular filter having a second diameter which is about twicethe first diameter in the middle region.
 18. An apparatus for measuringlung volume from a segmented lung image obtained from a scan whichincludes a plurality of slices of voxel data, the voxel data having agray level value and a volume associated therewith, the apparatuscomprising: a lung volume measuring unit configured to: generate amatrix of entries from the segmented lung image, the matrix having aplurality of columns and a plurality of rows, each of the plurality ofcolumns representing the gray level value and each of the plurality ofrows representing one of the plurality of slices in the scan, theentries corresponding to the number of times that the gray level occursin the corresponding slice; determine a number of voxels in thesegmented lung image from the matrix entries; and multiply the number ofvoxels in the segmented lung image by the volume of each voxel.
 19. Anapparatus for measuring volume of tissue in a segmented lung imageobtained from a scan including a plurality of slices of voxel data, thevoxel data having a gray level value and a volume, the apparatuscomprising: a lung tissue measuring unit configured to: generate amatrix of entries from the segmented lung image, the matrix having aplurality of columns and a plurality of rows, each of the plurality ofcolumns representing the gray level value and each of the plurality ofrows representing one of the plurality of slices in the scan, theentries corresponding to the number of times that the gray level occursin the corresponding slice; determine a sum of tissue of voxels in thesegmented lung image from the matrix entries; and multiply the sum oftissue of voxels in the segmented lung image by the volume of eachvoxel.
 20. An apparatus for measuring volume of tissue as defined byclaim 19, wherein the sum of tissue of voxels is calculated by summingthe product of each matrix entry multiplied by the corresponding graylevel value divided by a gray level value assigned for tissue.
 21. Anarticle of manufacture for generating a lung mask for segmenting a lungimage from voxel data containing the lung image, noise, solid componentsand surrounding background information, the article comprising: amachine readable medium containing one or more programs which whenexecuted implement the steps of: (a) applying a median filter and a meanfilter to the voxel data to reduce noise; (b) thresholding the noisereduced voxel data to identify solid components and other structures;(c) identifying the surrounding background in the noise reduced voxeldata; (d) deleting the surrounding background from the noise reducedvoxel data; (e) labeling connected components of the voxel data; (f)determining the largest connected components to select a lung regionfrom the voxel data, the lung region having a geometric form; and (g)morphological filtering the voxel data to refine the geometric form ofthe lung region.
 22. An article of manufacture for generating a lungmask as defined by claim 21, wherein the voxel data is associated with aplurality of slices through a patient's lung, the slices beginning atabout the patient's shoulders; and step (a) is performed for the first25 percent of the plurality of slices only whereby the computation timeis substantially reduced.
 23. An article of manufacture for generating alung mask as defined by claim 21, wherein the median filter has a sizeof 4×4.
 24. An article of manufacture for generating a lung mask asdefined by claim 21, wherein the mean filter has a size of 1×3.
 25. Anarticle of manufacture for generating a lung mask as defined by claim21, wherein the voxel data in step (b) is thresholded at a gray level ofabout
 500. 26. An article of manufacture for generating a lung maskdefined by claim 21, wherein the largest connected components areassociated with more than about 1 percent of the voxel data.
 27. Anarticle of manufacture for generating a lung mask as defined by claim21, wherein the voxel data is associated with a plurality of slices, theslices being divided into a first end region, a middle region, and asecond end region; and the morphological filtering being performed with2D circular filter having a first diameter in the first end region andthe second end region; and the morphological filtering being performedwith 2D circular filter having a second diameter which is about twicethe first diameter in the middle region.
 28. An article of manufacturefor measuring lung volume from a segmented lung image obtained from ascan which includes a plurality of slices of voxel data, the voxel datahaving a gray level value and a volume associated therewith, the articlecomprising: a machine readable medium containing one or more programswhich when executed implement the steps of: generating a matrix ofentries from the segmented lung image, the matrix having a plurality ofcolumns and a plurality of rows, each of the plurality of columnsrepresenting the gray level value and each of the plurality of rowsrepresenting one of the plurality of slices in the scan, the entriescorresponding to the number of times that the gray level occurs in thecorresponding slice; determining a number of voxels in the segmentedlung image from the matrix entries; and multiplying the number of voxelsin the segmented lung image by the volume of each voxel.
 29. An articleof manufacture for measuring volume of tissue in a segmented lung imageobtained from a scan including a plurality of slices of voxel data, thevoxel data having a gray level value and a volume, the articlecomprising: a machine readable medium containing one or more programswhich when executed implement the steps of: generating a matrix ofentries from the segmented lung image, the matrix having a plurality ofcolumns and a plurality of rows, each of the plurality of columnsrepresenting the gray level value and each of the plurality of rowsrepresenting one of the plurality of slices in the scan, the entriescorresponding to the number of times that the gray level occurs in thecorresponding slice; determining a sum of tissue of voxels in thesegmented lung image from the matrix entries; and multiplying the sum oftissue of voxels in the segmented lung image by the volume of eachvoxel.
 30. An article of manufacture for measuring volume of tissue asdefined by claim 29, wherein the sum of tissue of voxels is calculatedby summing the product of each matrix entry multiplied by thecorresponding gray level value divided by a gray level value assignedfor tissue.
 31. A method for finding the location, P′, and size, r′, ofa pulmonary nodule in a high-resolution computed tomography (CT) image,the method comprising: (a) selecting a set of initial processingparameters including an initial location, P₁, an initial size, r₁, andtarget value, T; (b) computing an initial new location P_(i) of thenodule with a locator template function; (c) incrementally increasingsize, r_(i); (d) computing a new location P_(i) of the nodule with thelocator template function; (e) computing a size metric with a sizingtemplate function; (f) repeating steps (c) through (e) until the sizemetric is less than the target value; (g) returning the location, P′,and the size, r′, from the previous iteration of steps (c) through (e).32. A method for finding the location, P, and size, r, of a pulmonarynodule in a high-resolution computed tomography image, the methodcomprising: (a) windowing the image to ignore bone structures; (b)selecting a locator template function and a sizing template function;(c) selecting a set of initial processing parameters including: aninitial location, P; size, r; and termination criteria; (d) performing asearch to determine a maximum response of the locator template function;(e) determining a response of the sizing template function and comparingthe response to the termination criteria; (f) incrementally increasingthe size, r, only if the termination criteria has not been satisfied andrepeating steps d and e; and (g) outputting the location, P, and size,r, of the nodule.
 33. A method for finding the location and size of apulmonary nodule as defined in claim 32, wherein the windowing comprisesclipping the image at intensities over about
 1000. 34. A method forfinding the location and size of a pulmonary nodule as defined in claim32, wherein the locator template function is a Gaussian templatefunction.
 35. A method for finding the location and size of a pulmonarynodule as defined in claim 32, wherein the locator template function isa Laplacian of the Gaussian template function.
 36. A method for findingthe location and size of a pulmonary nodule as defined in claim 32,wherein the locator template function is a difference of Gaussianstemplate function.
 37. A method for finding the location and size of apulmonary nodule as defined in claim 32, wherein the template functionhas at least four parameters corresponding to the x-location,y-location, z-location, and radius.
 38. A method for finding thelocation and size of a pulmonary nodule as defined in claim 32, whereinthe initial location, P, is calculated from the image.
 39. A method forfinding the location and size of a pulmonary nodule as defined in claim32, wherein the initial location, P, is specified by a user.
 40. Amethod for finding the location and size of a pulmonary nodule asdefined in claim 32, wherein the search is performed by a hill climbingmethod.
 41. A method for finding the location and size of a pulmonarynodule as defined in claim 32, wherein the search is performed byPowell's method.
 42. A pulmonary nodule finding apparatus for findingthe location, P′, and size, r′, of a pulmonary nodule in ahigh-resolution computed tomography (CT) image, the pulmonary nodulefinding apparatus comprising: a nodule finding unit configured to: (a)select a set of initial processing parameters including an initiallocation, P₁, an initial size, r₁, and target value, T; (b) compute aninitial new location P_(i) of the nodule with a locator templatefunction; (c) incrementally increase size, r_(i); (d) compute a newlocation P_(i) of the nodule with the locator template function; (e)compute a size metric with a sizing template function; (f) repeat steps(c) through (e) until the size metric is less than the target value; (g)return the location, P′, and the size, r′, from the previous iterationof steps (c) through (e).
 43. A pulmonary nodule finding apparatus forfinding the location, P, and size, r, of a pulmonary nodule in ahigh-resolution computed tomography image, the pulmonary nodule findingapparatus comprising: a nodule finding unit configured to: (a) windowthe image to ignore bone structures; (b) select a locator templatefunction and a sizing template function; (c) select a set of initialprocessing parameters including: an initial location, P; size, r; andtermination criteria; (d) perform a search to determine a maximumresponse of the locator template function; (e) determine a response ofthe sizing template function and comparing the response to thetermination criteria; (f) incrementally increase the size, r, only ifthe termination criteria has not been satisfied and repeating steps dand e; and (g) output the location, P, and size, r, of the nodule.
 44. Apulmonary nodule finding apparatus as defined in claim 43, wherein theimage is clipped at intensities over about 1000 to window the image. 45.A pulmonary nodule finding apparatus as defined in claim 43, wherein thelocator template function is a Gaussian template function.
 46. Apulmonary nodule finding apparatus as defined in claim 43, wherein thelocator template function is a Laplacian of the Gaussian templatefunction.
 47. A pulmonary nodule finding apparatus as defined in claim43, wherein the locator template function is a difference of Gaussianstemplate function.
 48. A pulmonary nodule finding apparatus as definedin claim 43, wherein the template function has at least four parameterscorresponding to the x-location, y-location, z-location, and radius. 49.A pulmonary nodule finding apparatus as defined in claim 43, wherein theinitial location, P, is calculated from the image.
 50. A pulmonarynodule finding apparatus as defined in claim 43, wherein the initiallocation, P, is specified by a user.
 51. A pulmonary nodule findingapparatus as defined in claim 43, wherein the search is performed by ahill climbing method.
 52. A pulmonary nodule finding apparatus asdefined in claim 43, wherein the search is performed by Powell's method.53. An article of manufacture for finding the location, P′, and size,r′, of a pulmonary nodule in a high-resolution computed tomography (CT)image, the article comprising: a machine readable medium containing oneor more programs which when executed implement the steps of: (a)selecting a set of initial processing parameters including an initiallocation, P₁, an initial size, r₁, and target value, T; (b) computing aninitial new location P_(i) of the nodule with a locator templatefunction; (c) incrementally increasing size, r_(i); (d) computing a newlocation P_(i) of the nodule with the locator template function; (e)computing a size metric with a sizing template function; (f) repeatingsteps (c) through (e) until the size metric is less than the targetvalue; (g) returning the location, P′, and the size, r′, from theprevious iteration of steps (c) through (e).
 54. An article ofmanufacture for finding the location, P, and size, r, of a pulmonarynodule in a high-resolution computed tomography image, the articlecomprising: a machine readable medium containing one or more programswhich when executed implement the steps of: (a) windowing the image toignore bone structures; (b) selecting a locator template function and asizing template function; (c) selecting a set of initial processingparameters including: an initial location, P; size, r; and terminationcriteria; (d) performing a search to determine a maximum response of thelocator template function; (e) determining a response of the sizingtemplate function and comparing the response to the terminationcriteria; (f) incrementally increasing the size, r, only if thetermination criteria has not been satisfied and repeating steps d and e;and (g) outputting the location, P, and size, r, of the nodule.
 55. Anarticle of manufacture for finding the location and size of a pulmonarynodule as defined in claim 54, wherein the windowing comprises clippingthe image at intensities over about
 1000. 56. An article of manufacturefor finding the location and size of a pulmonary nodule as defined inclaim 54, wherein the locator template function is a Gaussian templatefunction.
 57. An article of manufacture for finding the location andsize of a pulmonary nodule as defined in claim 54, wherein the locatortemplate function is a Laplacian of the Gaussian template function. 58.An article of manufacture for finding the location and size of apulmonary nodule as defined in claim 54, wherein the locator templatefunction is a difference of Gaussians template function.
 59. An articleof manufacture for finding the location and size of a pulmonary noduleas defined in claim 54, wherein the template function has at least fourparameters corresponding to the x-location, y-location, z-location, andradius.
 60. An article of manufacture for finding the location and sizeof a pulmonary nodule as defined in claim 54, wherein the initiallocation, P, is calculated from the image.
 61. An article of manufacturefor finding the location and size of a pulmonary nodule as defined inclaim 54, wherein the initial location, P, is specified by a user. 62.An article of manufacture for finding the location and size of apulmonary nodule as defined in claim 54, wherein the search is performedby a hill climbing method.
 63. An article of manufacture finding thelocation and size of a pulmonary nodule as defined in claim 54, whereinthe search is performed by Powell's method.
 64. A method for registering3-d images of a pulmonary nodule from a high-resolution computedtomography (CT) scans, the images being in a floating point pixel-formatassociated with a 6-dimensional parameter space and including a firstimage (im₁) obtained at time-1 and a second image (im₂) obtained attime-2, the method comprising the steps of: (a) calculating initialrigid-body transformation parameters for a rigid-body transformation onthe first image (im₁); (b) determining the optimum rigid-bodytransformation parameters by calculating a registration metric betweenthe second image (im₂) and the rigid-body transformation on the firstimage (im₁); and (c) generating a registered image from the optimumrigid-body transformation parameters.
 65. A method as defined in claim64, wherein step (a) is preceded by masking one of the images by settingpixels to a background value.
 66. A method as defined in claim 65,wherein the background value is about −1000.
 67. A method as defined inclaim 64, wherein the registration metric is minimized.
 68. A method asdefined in claim 64, wherein the registration metric is maximized.
 69. Amethod as defined in claim 64, wherein the registration metric iscalculated by transforming the first image (im₁) with the initialrigid-body transformation parameters to obtain a transformed first image(im_(1t)); calculating the registration metric as a correlation (C)between the transformed first image (im_(1t)) and the second image(im₂); and searching for the maximum correlation (C) in the6-dimensional parameter space.
 70. A method as defined in claim 64,wherein the registration metric is calculated by transforming the firstimage (im₁) with the initial rigid-body transformation parameters toobtain a transformed first image (im_(1t)); calculating the registrationmetric as a mean-squared-difference (MSD) between the transformed firstimage (im_(1t)) and the second image (im₂); and searching for theminimum mean-squared-difference (MSD) in the 6-dimensional parameterspace.
 71. A method as defined in claim 69, wherein the transforming ofthe first image (im₁) to obtain the transformed first image (im_(1t)) isa mapping of a point v in 3-d space to a point v′ in transformed spacedefined by: $\nu^{\prime} = {{R_{x}R_{y}R_{Z}\nu} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}$

wherein R_(x), R_(y), and R_(z) are rotation matrices defined as:$\begin{matrix}{R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {{\quad \quad^{-}}{\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}} \\{R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{{\quad \quad^{-}}{\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}} \\{R_{z} = \begin{bmatrix}{\cos \left( r_{z} \right)} & {{\quad \quad^{-}}{\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}_{\bullet}}\end{matrix}$


72. A method as defined in claim 70, wherein the transforming of thefirst image (im₁) to obtain the transformed first image (im_(1t)) is amapping of a point v in 3-d space to a point v′ in transformed spacedefined by: wherein R_(x), R_(y), and R_(z) are rotation matricesdefined as: $\begin{matrix}{R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {{\quad \quad^{-}}{\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}} \\{R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{{\quad \quad^{-}}{\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}} \\{R_{z} = \begin{bmatrix}{\cos \left( r_{z} \right)} & {{\quad \quad^{-}}{\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}_{\bullet}}\end{matrix}$


73. A method as defined in claim 64, wherein initial rigid-bodytransformation parameters include six parameters (tx,ty,tz,rx,ry,rz)respectively defined as translation in x, translation in y, translationin z, rotation about the x-axis, rotation about the y-axis, and rotationabout the z-axis; wherein the initial rotation parameters (rx,ry,rz) areall set to zero; and the initial translation parameters (tx,ty,tz,) areset so that the nodule in the first image (im₁) overlaps the nodule inthe second image (im₂) during the initial calculation of theregistration metric.
 74. A method as defined in claim 73, wherein theinitial translation parameters (tx,ty,tz,) are set to a differencebetween the center of the first image (im₁) and the center of the secondimage (im₂).
 75. A method as defined in claim 73, wherein the initialtranslation parameters (tx,ty,tz) are set to a difference between thecenter of mass of the first image (im₁) and the center of mass of thesecond image (im₂).
 76. A method as defined in claim 69, wherein thesearching is conducted by calculating the correlation (C) for everypossible set of rigid-body transformation parameters.
 77. A method asdefined in claim 70, wherein the searching is conducted by calculatingthe mean-squared-difference (MSD) for every possible set of rigid-bodytransformation parameters.
 78. A method as defined in claim 69, whereinthe searching is conducted by a Hill-Climbing search method.
 79. Amethod as defined in claim 70, wherein the searching is conducted by aHill-Climbing search method.
 80. A method as defined in claim 69,wherein the searching is conducted by Powell's method.
 81. A method asdefined in claim 70, wherein the searching is conducted by Powell'smethod.
 82. A registering apparatus for registering 3-d images of apulmonary nodule from a high-resolution computed tomography (CT) scans,the images being in a floating point pixel-format associated with a6-dimensional parameter space and including a first image (im₁) obtainedat time-1 and a second image (im₂) obtained at time-2, the registeringapparatus comprising: a registering unit configured to: (a) calculateinitial rigid-body transformation parameters for a rigid-bodytransformation on the first image (im₁); (b) determine the optimumrigid-body transformation parameters by calculating a registrationmetric between the second image (im₂) and the rigid-body transformationon the first image (im₁); and (c) generate a registered image from theoptimum rigid-body transformation parameters.
 83. A registeringapparatus as defined in claim 82, wherein one of the images is initiallymasked by setting pixels to a background value.
 84. A registeringapparatus as defined in claim 83, wherein the background value is about−1000.
 85. A registering apparatus as defined in claim 82, wherein theregistration metric is minimized.
 86. A registering apparatus as definedin claim 82, wherein the registration metric is maximized.
 87. Aregistering apparatus as defined in claim 82, wherein the registrationmetric is calculated by transforming the first image (im₁) with theinitial rigid-body transformation parameters to obtain a transformedfirst image (im_(1t)); calculating the registration metric as acorrelation (C) between the transformed first image (im_(1t)) and thesecond image (im₂); and searching for the maximum correlation (C) in the6-dimensional parameter space.
 88. A registering apparatus as defined inclaim 82, wherein the registration metric is calculated by transformingthe first image (im₁) with the initial rigid-body transformationparameters to obtain a transformed first image (im_(1t)); calculatingthe registration metric as a mean-squared-difference (MSD) between thetransformed first image (im_(1t)) and the second image (im₂); andsearching for the minimum mean-squared-difference (MSD) in the6-dimensional parameter space.
 89. A registering apparatus as defined inclaim 87, wherein the transforming of the first image (im₁) to obtainthe transformed first image (im_(1t)) is a mapping of a point v in 3-dspace to a point v′ in transformed space defined by:$\nu^{\prime} = {{R_{x}R_{y}R_{Z}\nu} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}$

wherein R_(x), R_(y), and R_(z) are rotation matrices defined as:$\begin{matrix}{R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {{\quad \quad^{-}}{\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}} \\{R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{{\quad \quad^{-}}{\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}} \\{R_{z} = \begin{bmatrix}{\cos \left( r_{z} \right)} & {{\quad \quad^{-}}{\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}_{\bullet}}\end{matrix}$


90. A registering apparatus as defined in claim 88, wherein thetransforming of the first image (im₁) to obtain the transformed firstimage (im_(1t)) is a mapping of a point v in 3-d space to a point v′ intransformed space defined by:$v^{\prime} = {{R_{x}R_{y}R_{z}v} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}$

wherein R_(x), R_(y), and R_(z) are rotation matrices defined as:$R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {- {\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}$ $R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{- {\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}$ $R_{z} = {\begin{bmatrix}{\cos \left( r_{z} \right)} & {- {\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}.}$


91. A registering apparatus as defined in claim 82, wherein initialrigid-body transformation parameters include six parameters(tx,ty,tz,rx,ry,rz) respectively defined as translation in x,translation in y, translation in z, rotation about the x-axis, rotationabout the y-axis, and rotation about the z-axis; wherein the initialrotation parameters (rx,ry,rz) are all set to zero; and the initialtranslation parameters (tx,ty,tz,) are set so that the nodule in thefirst image (im₁) overlaps the nodule in the second image (im₂) duringthe initial calculation of the registration metric.
 92. A registeringapparatus as defined in claim 91, wherein the initial translationparameters (tx,ty,tz,) are set to a difference between the center of thefirst image (im₁) and the center of the second image (im₂).
 93. Aregistering apparatus as defined in claim 91, wherein the initialtranslation parameters (tx,ty,tz) are set to a difference between thecenter of mass of the first image (im₁) and the center of mass of thesecond image (im₂).
 94. A registering apparatus as defined in claim 87,wherein the searching is conducted by calculating the correlation (C)for every possible set of rigid-body transformation parameters.
 95. Aregistering apparatus as defined in claim 88, wherein the searching isconducted by calculating the mean-squared-difference (MSD) for everypossible set of rigid-body transformation parameters.
 96. A registeringapparatus as defined in claim 87, wherein the searching is conducted bya Hill-Climbing search method.
 97. A registering apparatus as defined inclaim 88, wherein the searching is conducted by a Hill-Climbing searchmethod.
 98. A registering apparatus as defined in claim 87, wherein thesearching is conducted by Powell's method.
 99. A registering apparatusas defined in claim 88, wherein the searching is conducted by Powell'smethod.
 100. An article of manufacture for registering 3-d images of apulmonary nodule from a high-resolution computed tomography (CT) scans,the images being in a floating point pixel-format associated with a6-dimensional parameter space and including a first image (im₁) obtainedat time-1 and a second image (im₂) obtained at time-2, the articlecomprising: a machine readable medium containing one or more programswhich when executed implement the steps of: (a) calculating initialrigid-body transformation parameters for a rigid-body transformation onthe first image (im₁); (b) determining the optimum rigid-bodytransformation parameters by calculating a registration metric betweenthe second image (im₂) and the rigid-body transformation on the firstimage (im₁); and (c) generating a registered image from the optimumrigid-body transformation parameters.
 101. An article of manufacture forregistering 3-d images of a pulmonary nodule as defined in claim 100,wherein step (a) is preceded by masking one of the images by settingpixels to a background value.
 102. An article of manufacture forregistering 3-d images of a pulmonary nodule as defined in claim 101,wherein the background value is about −1000.
 103. An article ofmanufacture for registering 3-d images of a pulmonary nodule as definedin claim 100, wherein the registration metric is minimized.
 104. Anarticle of manufacture for registering 3-d images of a pulmonary noduleas defined in claim 100, wherein the registration metric is maximized.105. An article of manufacture for registering 3-d images of a pulmonarynodule as defined in claim 100, wherein the registration metric iscalculated by transforming the first image (im₁) with the initialrigid-body transformation parameters to obtain a transformed first image(im_(1t)); calculating the registration metric as a correlation (C)between the transformed first image (im_(1t)) and the second image(im₂); and searching for the maximum correlation (C) in the6-dimensional parameter space.
 106. An article of manufacture forregistering 3-d images of a pulmonary nodule as defined in claim 100,wherein the registration metric is calculated by transforming the firstimage (im₁) with the initial rigid-body transformation parameters toobtain a transformed first image (im_(1t)); calculating the registrationmetric as a mean-squared-difference (MSD) between the transformed firstimage (im_(1t)) and the second image (im₂); and searching for theminimum mean-squared-difference (MSD) in the 6-dimensional parameterspace.
 107. An article of manufacture for registering 3-d images of apulmonary nodule as defined in claim 105, wherein the transforming ofthe first image (im₁) to obtain the transformed first image (im_(1t)) isa mapping of a point v in 3-d space to a point v′ in transformed spacedefined by: $v^{\prime} = {{R_{x}R_{y}R_{z}v} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}$

wherein R_(x), R_(y) and R_(z) are rotation matrices defined as:$R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {- {\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}$ $R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{- {\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}$ $R_{z} = {\begin{bmatrix}{\cos \left( r_{z} \right)} & {- {\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}.}$


108. An article of manufacture for registering 3-d images of a pulmonarynodule as defined in claim 106, wherein the transforming of the firstimage (im₁) to obtain the transformed first image (im_(1t)) is a mappingof a point v in 3-d space to a point v′ in transformed space defined by:$v^{\prime} = {{R_{x}R_{y}R_{z}v} + \begin{bmatrix}t_{x} \\t_{y} \\t_{z}\end{bmatrix}}$

wherein R_(x), R_(y) and R_(z) are rotation matrices defined as:$R_{x} = \begin{bmatrix}1 & 0 & 0 \\0 & {\cos \left( r_{x} \right)} & {- {\sin \left( r_{x} \right)}} \\0 & {\sin \left( r_{x} \right)} & {\cos \left( r_{x} \right)}\end{bmatrix}$ $R_{y} = \begin{bmatrix}{\cos \left( r_{y} \right)} & 0 & {\sin \left( r_{y} \right)} \\0 & 1 & 0 \\{- {\sin \left( r_{y} \right)}} & 0 & {\cos \left( r_{y} \right)}\end{bmatrix}$ $R_{z} = {\begin{bmatrix}{\cos \left( r_{z} \right)} & {- {\sin \left( r_{z} \right)}} & 0 \\{\sin \left( r_{z} \right)} & {\cos \left( r_{z} \right)} & 0 \\0 & 0 & 1\end{bmatrix}.}$


109. An article of manufacture for registering 3-d images of a pulmonarynodule as defined in claim 100, wherein initial rigid-bodytransformation parameters include six parameters (tx,ty,tz,rx,ry,rz)respectively defined as translation in x, translation in y, translationin z, rotation about the x-axis, rotation about the y-axis, and rotationabout the z-axis; wherein the initial rotation parameters (rx,ry,rz) areall set to zero; and the initial translation parameters (tx,ty,tz,) areset so that the nodule in the first image (im₁) overlaps the nodule inthe second image (im₂) during the initial calculation of theregistration metric.
 110. An article of manufacture for registering 3-dimages of a pulmonary nodule as defined in claim 109, wherein theinitial translation parameters (tx,ty,tz,) are set to a differencebetween the center of the first image (im₁) and the center of the secondimage (im₂).
 111. An article of manufacture for registering 3-d imagesof a pulmonary nodule as defined in claim 109, wherein the initialtranslation parameters (tx,ty,tz) are set to a difference between thecenter of mass of the first image (im₁) and the center of mass of thesecond image (im₂).
 112. An article of manufacture for registering 3-dimages of a pulmonary nodule as defined in claim 105, wherein thesearching is conducted by calculating the correlation (C) for everypossible set of rigid-body transformation parameters.
 113. An article ofmanufacture for registering 3-d images of a pulmonary nodule as definedin claim 106, wherein the searching is conducted by calculating themean-squared-difference (MSD) for every possible set of rigid-bodytransformation parameters.
 114. An article of manufacture forregistering 3-d images of a pulmonary nodule as defined in claim 105,wherein the searching is conducted by a Hill-Climbing search method.115. An article of manufacture for registering 3-d images of a pulmonarynodule as defined in claim 106, wherein the searching is conducted by aHill-Climbing search method.
 116. An article of manufacture forregistering 3-d images of a pulmonary nodule as defined in claim 105,wherein the searching is conducted by Powell's method.
 117. An articleof manufacture for registering 3-d images of a pulmonary nodule asdefined in claim 106, wherein the searching is conducted by Powell'smethod.
 118. A method for removing extraneous matter from an image, theimage including a juxtapleural nodule, the method comprising the stepsof: (a) providing an initial location P′; (b) calculating a sphericalvolume centered at the initial location P′, the spherical volume fittinginside the image; (c) calculating a center of mass COM of the sphericalvolume; (d) determining an initial direction d′, the initial directiond′ being directed towards the extraneous matter in accordance with thefollowing equation${d^{\prime} = \frac{{COM} - P^{\prime}}{{{COM} - P^{\prime}}}};$

(e) initializing a current location P_(i) to be equal to the initiallocation P′; (f) initializing a current direction di to be equal to theinitial direction d′; (g) initializing a maximum ratio γ_(max), stepsize s, prior mass mass_(i-1), and prior change in mass Δ_(i-1); (h)moving the current location P_(i) by the step size s in the currentdirection d_(i); (i) determining an equation defining a plane A, theplane A being normal to the current direction d_(i), the plane A passingthrough the current location P_(i); (j) calculating a current massmass_(i) of the nodule on a side of the plane A opposing that of theextraneous matter; (k) calculating a current change in mass Δ_(i) bysubtracting the prior mass mass_(i-1) from the current mass mass_(i);(l) calculating a current ratio γ in accordance with the followingequation ${\gamma = {\frac{\Delta_{i}}{\Delta_{i - 1}} - 1}};$

(m) setting the prior mass mass_(i-1) equal to the current massmass_(i); (n) setting the prior change in mass Δ_(i-1) equal to thecurrent change in mass Δ_(i); (o) comparing the current ratio γ to themaximum ratio γ_(max); (p) modifying the current direction d_(i) tominimize the current mass mass_(i), and performing steps (h)-(o) whilethe current ratio γ is one of less than and equal to the maximum ratioγ_(max); and (q) modifying the current direction d_(i), performing steps(h)-(o), and outputting the area of the nodule partitioned by the planeA in response to the current ratio γ being greater than the maximumratio γ_(max).
 119. A method for removing extraneous matter from animage, the image including a juxtapleural nodule as defined in claim118, wherein the extraneous matter includes a pleural surface.
 120. Amethod for removing extraneous matter from an image, the image includinga juxtapleural nodule as defined in claim 118, wherein step (g) furtherincludes the steps of: initializing the maximum ratio γ_(max) to 0.5;and initializing the step size s to 1.5.
 121. A method for removingextraneous matter from an image, the image including a juxtapleuralnodule as defined in claim 118, wherein step (j) further includes thesteps of: (r) defining the current location P_(i) as being visited; (s)determining on which side of the plane A the current location P_(i) islocated; (t) terminating in reponse to the current location P_(i) notbeing located on a side of the plane opposing that of the extraneousmatter; (u) defining the current location P_(i) as being part of aregion of interest in response to the current location P_(i) beinglocated on the side of the plane opposing that of the extraneous matter;and (v) performing steps (r)-(u) recursively using a locationcorresponding to at least one of six (6) one-pixel moves from thecurrent location P_(i).
 122. A method for removing extraneous matterfrom an image, the image including a juxtapleural nodule as defined inclaim 118, wherein step (p) further includes the steps of: calculatingrecursively the current mass mass_(i) of the nodule on a side of theplane A opposing that of the extraneous matter using at least one of six(6) directions and a step size s₁ from the current location P_(i); anddefining the current direction d_(i) equal to the direction yielding thelargest decrease in the current mass mass_(i).
 123. A method forremoving extraneous matter from an image, the image including ajuxtapleural nodule as defined in claim 118, wherein step (q) furtherincludes the steps of: calculating recursively the current mass mass_(i)of the nodule on a side of the plane A opposing that of the extraneousmatter using at least one of six (6) directions and a step size s₁ fromthe current location P_(i); and defining the current direction d_(i)equal to the direction yielding the largest decrease in the current massmass_(i).
 124. A method for removing extraneous matter from an image,the image including a juxtapleural nodule as defined in claim 118,wherein the initial location P′ is located near a center of the nodule.125. A recursive apparatus for removing extraneous matter from an image,the image including a juxtapleural nodule, the recursive apparatuscomprising: a processing unit, the processing unit being configured to:(a) accept an initial location P′; (b) calculate a spherical volumecentered at the initial location P′, the spherical volume fitting insidethe image; (c) calculate a center of mass COM of the spherical volume;(d) determine an initial direction d′, the initial direction d′ beingdirected towards the extraneous matter in accordance with the followingequation${d^{\prime} = \frac{{COM} - P^{\prime}}{{{COM} - P^{\prime}}}};$

(e) initialize a current location P_(i) to be equal to the initiallocation P′; (f) initialize a current direction d_(i) to be equal to theinitial direction d′; (g) initialize a maximum ratio γ_(max), step sizes, prior mass mass_(i-1), and prior change in mass Δ_(i-1); (h) move thecurrent location P_(i) by the step size s in the current directiond_(i); (i) determine an equation defining a plane A, the plane A beingnormal to the current direction d_(i), the plane A passing through thecurrent location P_(i); (j) calculate a current mass mass_(i) of thenodule on a side of the plane A opposing that of the extraneous matter;(k) calculate a current change in mass Δ_(i) by subtracting the priormass mass_(i-1) from the current mass mass_(i); (l) calculate a currentratio γ in accordance with the following equation${\gamma = {\frac{\Delta_{i}}{\Delta_{i - 1}} - 1}};$

(m) set the prior mass mass_(i-1) equal to the current mass mass_(i);(n) set the prior change in mass Δ_(i-1) equal to the current change inmass Δ_(i); (o) compare the current ratio γ to the maximum ratioγ_(max); (p) modify the current direction d_(i) to minimize the currentmass mass_(i), and performing steps (h)-(o) while the current ratio γ isone of less than and equal to the maximum ratio γ_(max); and (q) modifythe current direction d_(i), performing steps (h)-(o), and outputtingthe area of the nodule partitioned by the plane A in response to thecurrent ratio γ being greater than the maximum ratio γ_(max).
 126. Arecursive apparatus for removing extraneous matter from an image asdefined in claim 125, wherein the extraneous matter includes a pleuralsurface.
 127. A recursive apparatus for removing extraneous matter froman image as defined in claim 125, wherein in step (g) the processingunit is configured to: initialize the maximum ratio γ_(max) to 0.5; andinitialize the step size s to 1.5.
 128. A recursive apparatus forremoving extraneous matter from an image as defined in claim 125,wherein in step (j) the processing unit is configured to: (r) define thecurrent location P_(i) as being visited; (s) determine on which side ofthe plane A the current location P_(i) is located; (t) terminate inreponse to the current location P_(i) not being located on a side of theplane opposing that of the extraneous matter; (u) define the currentlocation P_(i) as being part of a region of interest in response to thecurrent location P_(i) being located on the side of the plane opposingthat of the extraneous matter; and (v) perform steps (r)-(u) recursivelyusing a location corresponding to at least one of six (6) one-pixelmoves from the current location P_(i).
 129. A recursive apparatus forremoving extraneous matter from an image as defined in claim 125,wherein in step (p) the processing unit is configured to: calculaterecursively the current mass mass_(i) of the nodule on a side of theplane A opposing that of the extraneous matter using at least one of six(6) directions and a step size s from the current location P_(t1); anddefine the current direction di equal to the direction yielding thelargest decrease in the current mass mass_(i).
 130. A recursiveapparatus for removing extraneous matter from an image as defined inclaim 125, wherein in step (q) the processing unit is configured to:calculate recursively the current mass mass_(i) of the nodule on a sideof the plane A opposing that of the extraneous matter using at least oneof six (6) directions and a step size s₁ from the current locationP_(i); and defining the current direction d_(i) equal to the directionyielding the largest decrease in the current mass mass_(i).
 131. Arecursive apparatus for removing extraneous matter from an image asdefined in claim 125, wherein the initial location P′ is located near acenter of the nodule.
 132. An article of manufacture for removingextraneous matter from an image, the image including a juxtapleuralnodule, the article comprising: a machine readable medium containing atleast one program which when executed implements the steps of: (a)accepting an initial location P′; (b) calculating a spherical volumecentered at the initial location P′, the spherical volume fitting insidethe image; (c) calculating a center of mass COM of the spherical volume;(d) determining an initial direction d′, the initial direction d′ beingdirected towards the extraneous matter in accordance with the followingequation${d^{\prime} = \frac{{COM} - P^{\prime}}{\left. ||{{COM} - P^{\prime}} \right.||}};$

(e) initializing a current location P_(i) to be equal to the initiallocation P′; (f) initializing a current direction d_(i) to be equal tothe initial direction d′; (g) initializing a maximum ratio γ_(max), stepsize s, prior mass mass_(i-1), and prior change in mass Δ_(i-1); (h)moving the current location P_(i) by the step size s in the currentdirection d_(i); (i) determining an equation defining a plane A, theplane A being normal to the current direction d_(i), the plane A passingthrough the current location P_(i); (j) calculating a current massmass_(i) of the nodule on a side of the plane A opposing that of theextraneous matter; (k) calculating a current change in mass Δ_(i) bysubtracting the prior mass mass_(i-1) from the current mass mass_(i);(l) calculating a current ratio γ in accordance with the followingequation ${\gamma = {\frac{\Delta_{i}}{\Delta_{i - 1}} - 1}};$

(m) setting the prior mass mass_(i-1) equal to the current massmass_(i); (n) setting the prior change in mass Δ_(i-1) equal to thecurrent change in mass Δ_(i); o) comparing the current ratio γ to themaximum ratio γ_(max); (p) modifying the current direction d_(i) tominimize the current mass mass_(i), and performing steps (h)-(o) whilethe current ratio γ is one of less than and equal to the maximum ratioγ_(max); and (q) modifying the current direction d_(i), performing steps(h)-(o), and outputting the area of the nodule partitioned by the planeA in response to the current ratio γ being greater than the maximumratio γ_(max).
 133. An article of manufacture for removing extraneousmatter from an image as defined in claim 132, wherein the extraneousmatter includes a pleural surface.
 134. An article of manufacture forremoving extraneous matter from an image as defined in claim 132,wherein in step (g) the article further implements the steps of:initializing the maximum ratio γ_(max) to 0.5; and initializing the stepsize s to 1.5.
 135. An article of manufacture for removing extraneousmatter from an image as defined in claim 132, wherein in step (j) thearticle further implements the steps of: (r) defining the currentlocation P_(i) as being visited; (s) determining on which side of theplane A the current location P_(i) is located; (t) terminating inreponse to the current location P_(i) not being located on a side of theplane opposing that of the extraneous matter; (u) defining the currentlocation P_(i) as being part of a region of interest in response to thecurrent location P_(i) being located on the side of the plane opposingthat of the extraneous matter; and (v) performing steps (r)-(u)recursively using a location corresponding to at least one of six (6)one-pixel moves from the current location P_(i).
 136. An article ofmanufacture for removing extraneous matter from an image as defined inclaim 132, wherein in step (p) the article further implements the stepsof: calculating recursively the current mass mass_(i) of the nodule on aside of the plane A opposing that of the extraneous matter using atleast one of six (6) directions and a step size s₁ from the currentlocation P_(i); and defining the current direction d_(i) equal to thedirection yielding the largest decrease in the current mass mass_(i).137. An article of manufacture for removing extraneous matter from animage as defined in claim 132, wherein in step (q) the article furtherimplements the steps of: calculating recursively the current massmass_(i) of the nodule on a side of the plane A opposing that of theextraneous matter using at least one of six (6) directions and a stepsize s₁ from the current location P_(i); and defining the currentdirection d_(i) equal to the direction yielding the largest decrease inthe current mass mass_(i).
 138. An article of manufacture for removingextraneous matter from an image as defined in claim 132, wherein theinitial location P′ is located near a center of the nodule.